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SUMMARY 


This report describes a digital computer code CAVE (Conduction Analysis Via 
Eigenvalues), which finds application in the analysis of two-dimensional transient 
heating of hypersonic vehicles. The code is an extension of the work reported in NASA 
CR-2435 for the inverse conduction problem . CAVE is written in FORTRAN IV and 
is operational on both IBM 360-67 and CDC 6600 computers. 

The main advantages of CAVE over more conventional thermal analyzer codes 

are: 

• The method of solution is a hybrid analytical-numerical technique that is 
inherently stable permitting large tims steps even with the best of 
conductors having the finest of mesh size. This method can provide a 
facto r-of-five reduction in machine time compared to conventional 
explicit finite difference methods when structures with small time constants 
are analyzed over long flight trajectories. 

• The aerodynamic heating boundary conditions are calculated by the code 
based on the input flight trajectory (i.e. , altitude, velocily and angle of 
attack as functions of time) rather than calculated external to the code 
and then entered as input data. 

• The code computes the network conduction and convection links, as well 
as capacitance values, given basic geometrical and mesh sizes, for four 
geometries (leading edges , cooled panels , X-24C structure and slabs) . 

• The output from the code at each time interval includes the steady-state 
solution corresponding to the boundary conditions for that time interval. 

• The code also permits direct input of the heat transfer couplings, node 
capacitances and boundary conditions . 


This report is primarily a user's manual for the CAVE code. Input and output 
formats are presented and explained. Sample problems are included. A brief 
summary of the hybrid analytical -numerical technique, which utilizes eigenvalues 
(thermal frequencies) and eigenvectors (thermal mode vectors) is given in an 
appendix. Other appendixes include the aerod5mamic heating equations that have been 
incorporated in the code and flow charts . 
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Section 1 


INTRODUCTION 

The computer code CAVE (Conductive Analysis Via Eigenvalues) provides a 
very convenient and efficient tool for predicting the temperatures within thermal 
protection systems for hypersonic vehicles. 

The CAVE code is convenient to use because, first, the boundary conditions 
(convective heat transfer coefficient and adiabatic wall temperature) are calculated 
by the code based on the input values for altitude, velocity and angle of attack as 
functions of time. And, second, because the conduction and convection links between 
nodes, and the capacitance of each node are calculated by the code for leading edges, 
slabs and cooled panels (Fig. 1). The code also permits direct input of the heat 
transfer couplings, node capacitances and boundary conditions for other two- 
dimensional problems. 

CAVE can be very efficient in the use of computer time because the method 
employed to solve the partial differential heat conduction equation is a hybrid 
analytical-numerical (HAN) technique. In this method, spatial derivatives are 
replaced by appropriate finite difference representations and the temporal deriva- 
tives are retained as ordinary derivatives. In effect the problem is subdivided into 
a number of uniform temperature systems or nodes that are coupled and changing in 
temperature. The problem is thereby specified by a set of first order, linear, 
ordinary differential equations. The solution to the set of equations is expressed in 
terms of eigenvectors (thermal mode vectors for the system) and eigenvalues (thermal 
frequencies of the system). Appendix A gives details of the method. The important 
thing to note is that this method is particularly efficient in the use of computer time 
when the heat flux response is contained in the first few thermal modes (character- 
istic of materials with high thermal diffusivity) or if the response for a large number 
of time increments is required, which is precisely the situation in predicting the 
temperatures throughout the flight trajectory of a hypersonic vehicle. A reduction by 
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FIG. 1 GEOMETRIES BUILT INTO CAVE CODE 
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a factor of five in computer time can be expected over conventional explicit finite 
difference codes for typical flight trajectory analyses. The savings in computer 
time is due to the HAN method being inherently stable and, therefore, permitting 
large time steps,. 

This report is basically a user's manual for CAVE. Section 2 describes the 
overall operation and running of CAVE while Section 3 discusses the details of the 
leading edge problem. The input data format is presented and the output from sample 
problems is reviewed. Section 4 provides a similar treatment for the slab, cooled 
panel and X-24C geometries. Section 5 discusses an arbitrarily shaped geometry. 
The appendixes present such information as the details of the HAN method (Appendix 
A), the aerodynamic heating equations (Appendix B), a discussion of the treatment 
of radiation (Appendix C), a brief description of the CAVE subroutines (Appendix D), 
a discussion of nonlinearities and time dependency of h and T^^ (Appendix E), and 
a derivation of the eigenvalue/eigenvector solution (Appendix F). 

Mr, James L. Hunt, of the High Speed Aerodynamics Division, Langley Re- 
search Center, Virginia, served as the NASA technical monitor for the program. 

At Grumman, the contract was administered by the Advanced Development 
office, imder Mr. Fred Berger, Manager of Advanced Development System Engin- 
eering, The Study Manager was Dr. Kenneth A, Rathjen. 

Mr. Michael J. Rossi served as numerical analysis consultant for the pro- 
gram. Mr. Rossi developed the numerical method and the matrix subroutines pack- 
age under contract NAS 1-11818, "Lateral Conduction Effects on Heat- Transfer Data 
Obtained with the Phase- Change Paint Technique, " described in report NASA 
CR-2435 with co-author George Maise. Mr. Hunt also served as technical monitor 
for that program. 

Messrs. William Timlen and Charles Osonitsch were of considerable assist- 
ance in providing the appropriate aerodjmamic heating functions. Mr. Timlen also 
gave important support by running Grumman' s TTAl computer code to obtain in- 
dependent checks on the CAVE code. Mr. Brian Martin developed subroutine 
X-24C. 

The many helpful discussions with Dr. Gianky DaFomo are gratefully ac- 
knowledged. 


5 




Section 2 


DESCRIPTION OF CAVE CODE OPERATION 

This section provides an overview of the CAVE code capabilities, input/ output, 
and method of solution. 

CAVE has been designed with the convenience of the user in mind. Usual opera- 
tion of CAVE requires only the following elemental information from the user: 

• Selection of one of the built-in configurations of Figure 1 or the general 
geometry option 

• Specification of the surface emissivity and background radiation tempera- 
ture for problems involving radiation heat transfer 

• Specification of the material density, specific heat and thermal conductivity 
(the latter two of which can be temperature dependent) 

• Geometry-type information such as overall dimensions and grid network 
sizes 

• Initial temperature distribution 

• Flight trajectory, i.e. , altitude, velocity and angle of attack as functions 

of time in tabular form. This information is used by the code to predict 
the aerodynamic heat transfer coefficient, h, and adiabatic wall tempera- 
ture, as functions of time. Optionally the user may supply tables for 

h and T as functions of time and distance 

AW 

• Specification of the time step intervals 

With the above information specified, the problem solution is accomplished in 
the following sequence by the code: 

1. Storage requirements for the various arrays are determined and allocated. 

2. The geometry is discretized into elements and the volumes, conduction areas 
and lengths are computed. (A unit depth is assumed by the code.) For the 
general geometry problem these quantities are input data. 
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Then for each time step the following are done; 

3. Using the temperature distribution at the start of the time step, the 
thermal properties of the materials are determined, followed by the 
capacitances and conductances for the network. This step is exercised 
just once if the material properties are independent of temperature. 

4. Using the flight trajectory data at the start and end of the time step, the 
code establishes the time average convective heat transfer coefficient 
at each surface node and the corresponding time average adiabatic wall 
temperatures. Appendix B gives the particular aerodynamic heating 
equations used. 

5. The convective heat transfer couplings, due either to aerodynamic heating 
or internal coolant flow, at each surface node are then determined by taking 
the product of the convective coefficient and the surface area. These coup- 
lings are then modified to account for radiation if it is being considered . 
Appendix C gives the details of the linearized treatment that is given to the 
radiation heat transfer. 


At this jimcture we may visualize the code as being faced with the task of finding 
the solution to the following system of n first-order linear differential equations with 
constant coefficients ; 


C. i = K.. (T. - T.) + H. (T.„, . - T.) 

1 IT" U 3 1 1 AW, 1 


where = thermal capacitance of node i 

= conductive coupling between nodes i and j 
T^ - temperature of node i 

T^ = temperature of node j which is adjacent to node i 

H. = convective coupling between node i and the fluid (for interior 
nodes = 0) 

t == time 

T . = adiabatic wall temperature of the fluid in contact with node i 
AW,i 

There are n such coupled differential equations, one for each of the n nodes. 
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It may be Interesting to digress for a moment to note that the usual thermal 

dT. 

analyzers take equation (1) a step further and replace the ordinary derivative ^ 

with a finite difference approximation. Depending on the form of the approximation 
either an explicit or implicit algorithm is obtained. In the common explicit and 
implicit formulations, the T^^’s and T.'s are taken to be constant during the time step 
Interval. In the current HAN method, the ordinary derivative is retained and the T^^'s 
and Tj's are treated as time -dependent variables in Eq. (1). This leads to a more 
accurate solution with no limitation on the time step from a stability standpoint. 
However, in solving Eq. (1), the HAN method treats the C., K.., H. and T . as 
constants. This is necessary, as discussed in Appendix E, for an eigenvalue 
solution to exist. The technique used within CAVE to handle variations in these 
parameters is to subdivide the total time interval (i.e. , the flight trajectory) and 
take these parameters to be piecewise constant within each time subinterval. 

Thus, the single problem of determining the temperature distribution in the 
structure for the entire flight trajectory where the boundary conditions are varying is 
solved by considering a number of subproblems where the boundary conditions are 
piecewise constant. These subproblems are interconnected in that the temperature at 
the end of one time subinterval becomes the initial temperature for the next time 
subinterval. It should be noted that the time subintervals, or time steps in the HAN 
method, are ts^ically of the order of seconds or tens of seconds which is probably 
100 to 1000 times larger than is permissible with the explicit method. 

CAVE arithmetically averages the convective coefficient and adiabatic tem- 
perature at the beginning and end of the time interval. Therefore in selecting the time 
subintervals, the user should be guided by the variation in the flight trajectory with 
particular concern for abrupt changes that affect the convective heating. Assuming 
the flight trajectory table has been set up with these important points of change, time 
subintervals equal to those used in the flight trajectory table will very often prove 
satisfactory. For those problems in which the temperature dependency of the 
material properties plays a dominant role for some reason , or if radiation heat 
transfer is of great importance, a second run with smaller subintervals should 
be made to determine the effect of sub interval selection on the predicted tem- 
peratures. 


The system of equations given in Eq. (1) has the following exact solution for a 


particular time -subinterval (refer to Appendix F): 


T. 

1 



c.. exp (X. t) 


Eq. (2) 


where 


T. 

1 


T 


00 


i 


c.. 


1 ] 


t 


= temperature at node i at time t seconds into the time subinterval 
= steady-state temperature at node i for the particular time subinterval 

= constants that depend on the , a set of eigenvectors of a matrix A, 

i 

and the temperatures of the nodes at the start of the time subinterval 
= the eigenvalues of a matrix A 

= time into the particular time subinterval. If T represents the time in the 
flight trajectory, and if t and T represent the time at the start and 
end of a time subinterval, then the followli^ relationships hold; 

0<t<T - T and T = T + t for T ^t<T 
es s s"“e 


A 


symmetric matrix whose elements depend on the C^, and H. of 
Eq. (1) , (Refer to Appendixes A and F) 


Considering a thermal network with 100 nodes, there are then 100 eigenvalues and 
eigenvectors to be determined and used in Eq. (2) . Considerable machine time can 
be saved by calculating only those eigenvalues and eigenvectors that are "significant” 
or "dominant". This was noted very aptly by Maise and Rossi in NASA CR-2435 and 
used by them in the CAPE code for the inverse heat transfer problem of finding the 
boundary conditions given the temperature history. When the series in Eq. (2) is 
truncated to the "dominant" terms, we obtain: 


T. 

1 


ne 



c.. exp (Xjt) 


Eq. (3) 
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where ne is a number substantially less than n. It represents the number of dominant 
eigenvalues and eigenvectors that will be foimd and utilized by the code. This is an 
input number decided upon by the user. Values of ne from 3 to 5 are recommended 
for most problems. Appendix A presents a discussion on this subject. 

With this background information, we are in a position to continue with the 
sequence that the CAVE code undergoes in finding the temperature history of the struc- 
ture throughout the flight trajectory. The next four steps involve matrix subroutines 
which for the most part were written and developed under contract NAS 1-11818 by 
M. J, Rossi. The sequence then continues from p. 6 as follows: 

6. Set up a matrix A in compact form which depends on the K.^ and H. of 
Eq. (1) (Refer to Appendix F). 

7. Obtain the ne dominant eigenvectors and eigenvalues of matrix A using 
Jennings* method of simultaneous vector iteration. 

8. Determine the steady-state solution to Eq. (1). 

9. Calculate the c^^ of Eq. (3) for i = 1, 2, . . . ,n and j = 1, 2, . . . ,ne. 

10. Calculate the temperatures of the nodes at the end of the time subinterval 
using Eq. (3). 

11. Set the initial temperatures for the next time subinterval equal to the final 
temperature of the present subinterval. Increment the flight trajectory 
time. 

12. Repeat steps 3 through 11 until the final time has been reached. The 
solution is then completed. 


The output from CAVE is for the most part self-explanatory and ’will be reviewed 
in detail in the following sections when sample problems are considered. In essence, 
there are three sections to the output. First, there is a partial feedback of the input 
data, including: geometric parameters, material properties, flight trajectory, or 
convective heat -transfer coefficient and adiabatic wall temperature, as the case may 
be, and initial temperature distribution. Secondly, there are the node numbers, 
material numbers, capacitances and conductances that were calculated by the code 


*A, Jennings, "A Direct Iteration Method of Obtaining Latent Roots and Vectors of a 
Symmetric Matrix," Proc. Cambridge Phil. Soc., 63, 1967, pp. 755-765, 
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for each node. And thirdly, there are printed out for each time subinterval the time 
average heat transfer coefficients , heat tr^sfer couplings and adiabatic wall temper- 
atures. Also printed out are Mach number, altitude, velocity, angle of attack and 
node temperatures at the end of each subinterval. As an aid to the user for a better 
feel of the problem being analyzed, the following items are also printed out: the 
steady-state temperature distribution for the time subinterval and the time -integrated 
heat input to each boundary node. 

Although CAVE has been designed to be most convenient for users interested 
in predicting structural temperatures of hypersonic vehicles, the code also proves 
convenient for analyzing the geometries given in Figure 1, subjected to other than the 
normal aerodynamic heating. The user may take advantage of the automatic division 
of these geometries by the code and supply the particular boundary conditions of his 
problem as input data. CAVE also proves a valuable code for analyzing geometries 
other than those given in Figure 1, i.e. , whenever the time constant of the system is 
small compared to total time of interest. In this case, the HAN method of CAVE 
offers significant machine-time savings compared to conventional methods. Sections 
3 and 4 consider in detail the built-in geometries and contain sample cases. Section 5 
discusses the general two-dimensional capabilities of CAVE. 
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Section 3 


LEADING EDGE GEOMETRY 


3.1 DISCUSSION 

This section presents the leading edge geometry that has been incorporated into 
CAVE and discusses how this geometry is discretized into nodes by the LEAD4 sub- 
routine. * This section also presents the input data format for this geometry. 

Figure 2 shows the leading edge geometry that is tacitly assumed by subroutine 
LEAD4 when it generates a nodal network. The insulating layer can be eliminated, 
as can either one or both of the coolant passages. These eliminations are accom- 
plished very simply by using input values of zero for the insulation thickness and 
coolant passage radii. 

Figure 3 shows the grid network for the leading edge. The nose region is 
divided into elements by concentric arcs and rays . The wedge portion is divided into 
rectangular elements, except near the coolant passages where odd shapes are 
encountered, and near the centerline, where the elements are trapezoidal. The cal- 
culation of the node capacitances is a straight-forward matter and it is done exactly 
for all elements including the trapezoids. The calculation of the conductances is also 
a straight-forward matter except near the coolant passages where the conductances 
are approximated using an effective area and length between nodes. In regard to 
the approximation, the code requires that the following relationships be maintained 
(Fig. 3): 

AXf = RPl, radius of nose coolant passage 

= ^^i+2 " T where RP2 = radius of aft coolant passage 

(and so, for example, i = 4 in Fig. 3) 

*LEAD4 is an expanded version of LEAD, which is a subroutine for leading edges 
written by George Maise under contract NAS 1-11818 and reported in NASA CR-2435 
by George Maise and Michael J. Rossi. The expanded version can handle cooling 
passages and a layer of insulation applied around the leading edge. 
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It is not necessary that the three AX's associated with the aft coolant passage 
be numbers 4, 5, and 6; that is just how it worked out in Figure 3. If there is no 
nose coolant passage, then Ax^ can be arbitrary in size (but nonzero). 

The thickness of the insulating material is given by TAU and it may equal zero. 
Notice in Figure 3 that nodes are located at this interface between the two. materials. 
LEAD4 assumes that there are equal volumes of the two materials associated with 
each interface node. Meaning that one half of Ar^, in this case, is associated with 
the insulator and the other half with the main material . 

The user may elect to have CAVE calculate the convective heat transfer coef- 
ficients and adiabatic wall temperatures around the leading edge, or he may supply 
tabular inputs for them. If the user elects the former option, then he supplies tabular 
values for the flight parameters of velocity, altitude and angle of attack as functions 
of time; moreover, he flags CAVE to use either the turbulent or laminar flow 
correlations, the details of which are presented in Appendix B. 

For leading-edge problems that involve increased heating due to local interfer- 
ence heating or some other effect such as plume impingement during a portion of the 
flight trajectory, the user may input two tables into CAVE. The tabular values are 
multiplicative factors which are position and time dependent. Values from one table 
are used to modify the convective coefficient on the top surface and values from the 
other table are for the bottom surface. A nonzero value for the input variable 
HMODI flags CAVE that this heating multiplier option will be exercised. For the 
normal run when the convective coefficient is not to be modified, HMODI equals 0, 
and the tables are omitted. 

In using this multiplier option, it is important to bear in mind that in modifying 
the convective coefficients, CAVE takes the average of the multiplicative factor at 
the beginning and end of the time interval and applies it over the entire interval. 
Therefore, the tables and computing time intervals must be selected with same care 
whenever step changes are to be simulated. A sample problem in Section 4 illustrates 
this. 
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The node numbering scheme for the leading edge geometry is interesting, A 
review of Figure 3 shows the following pattern: 




When we get to the sample problem, we shall see that the temperatures and 
other nodal properties are printed out in the following array form: 


!• 
2% 

4 « 
5« 
6 « 
?• 
8 » 
9e 
K lOm 



13« 

25« 

14» 

26» 

15« 

etc. 

16« 


17« 


18* 


19» 


20» 


21* 


22« 


23» 


24* 
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Two observations can be made. First, the nodes along the top surface of the 
leading edge (numbers 1, 13, 25, . . .) are printed out as the first row of the array 
and the nodes along the bottom surface (numbers 12, 24, 36, . . .) are printed out 
as the last row. And, secondly, if the elements of the first column are rotated as the 
arrows indicate, the array gets rearranged into something looking somewhat like the 
nodal arrangement within the leading edge. With a little experience, the user of 
CAVE is able to quickly scan the output and get an immediate feel for the temperature 
gradients within the leading edge. 

The following subsections present the input data format for leading-edge problems 
as well as for a sample problem . 

3.2 INPUT DATA FORMAT FOR LEADING EDGE GEOMETRY 
Indexes Card 

• JGEO, L, M, NE* (415) 

JGEO = 1 (selects leading edge geometry) 

L = number of elements through the material (must be an 
even integer) 

M = number of elements along top (or bottom) half of 
leading edge 

NE ~ number of dominant eigenvalues to be used in solution 
(e.g. , a typical number is 5) 

Title Card 

• Run identification, comments, etc. (5A10) 

Radiation Card 

• EPSl, TBGl (2F10.5) 

EPSl = emissivity of surface 

TBGl = background radiation temperature, °R 


♦The product L x M equals the total number of nodes. The current dimension 
statements in CAVE require that L x M<200 and that M < 25 for the leading edge 
geometry. 
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Material Properties Cards 


• 

MAT 



(15) 

• 

NMATl, RHOl, CONAVl, CPAVl 


(I10,3F10.5) 

• 

TMATl, TMAT1(2), ... 

, TMATl (NMATl) i 

j 

omit 

(8E10.0) 

• 

CONDTl(l), CONDTl(2), 

CONDTl (NMATl) * 

if 

(8E10.0) 

• 

CPTl(l), CPT1(2), ...» 

CPTl (NMATl) 

NMATl 

=0 (8E10.0) 

(If MAT 

= 2 include the cards:) 




• 

NMAT2, RH02, CONAV2, CPAV2 


(I10.0,3F10.5) 

• 

TMAT2(1), TMAT2(2), . 

. . , TMAT2 (NMAT2) ] 

omit 

(8E10.0) 

• 

CONDT2(l), CONDT2(2), 

, . . . , CONDT2 (NMAT2) 

if 

(8E10.0) 

• 

CPT2(1), CPT2(2), ..., 

CPT2 (NMAT2) ] 

NMAT2 

=0 (8E10.0) 


MAT 

NMATl 

RHOl 

CONAVl 

CPAVl 

TMATl (I) 


CONDTl (I) 
CPTl (I) 


= number of materials (1 or 2) 

= number of entries in properties table (maximum 
of 10). NMATl = 0 for constant properties 
= density of material 1, Ibm/cu-ft 
= average thermal conductivity of material 1 (used 
when NMATl = 0), Btu/ft-sec-®R 
= average specific heat of material 1 (used when 
NMATl = 0), Btu/lbm-°R 
= temperatures in thermal properties table for 
which CONDTl (I) and CPTl (I) are given; 

1 = 1, 2, . . NMATl, °R 
= thermal conductivity of material 1 at temperature 
TMATl (I), Btu/ft-sec-“R 
= specific heat of material 1 at temperatur6 
TMATl (I), Btu/lbm-“R 


NMAT2, RH02, CONAV2, etc. , same as NMATl, RHOl, CONAVl, etc. , 
except applied to material 2 
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Detail Geometry Cards 

• MCAP, THETA (I10.F10.5) 

• DELX(l), DELX(2), DELX(3), . . . , DELX(MM) (8F10.5) 

• DELR(l), DELR(2), DELR(3) DELR(L/2) (8F10.5) 

• TAH (F10.5) 

• RPl, RP2, S, HCOOLl, HCOOL2, TCOOLl, TCOOL2 (7F10.5) 


MCAP - number of elements into which nose of leading edge 
is subdivided (must be an even integer) 

THETA = wedge half angle of leading edge, in degrees 

DELX(I) = spatial increments in x direction 1=1, 2 MM 

(where MM = M - MC AP/2) , ft 

DELR(I) = spatial increments in radial direction 1=1, 2, , , , , L/2 
TAU = thickness of material 1, ft (when considering only one material, 
TAU = 0) 

RPl = radius of nose coolant passage, ft 
RP2 = radius of aft coolant passage, ft 
S = distance between coolant passage centers, ft 

HCOOLl = convective heat transfer coefficient inside nose coolant 
passage, Btu/ft -sec-TEl 

HCOOL2 = convective heat transfer coefficient inside aft coolant 
passage, Btu/ft -sec-®R 
TCOOLl = nose coolant temperature, °R 
TCOOL2 = aft coolant temperature, “R 

Initial Temperature Cards 

• KODE, I, T(I), H, JJ (215, ElO. 0,215) 

• 11100 (indicates end of initial temperature cards) (15) 

KODE = 0 or blank 
I = node number 
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1 

I 

I T(I) = node initial temperature, °R ' ■ 

I n and JJ = the node number is incremented by the spacing JJ until 

the limit II is reached. Each node so specified is 
c| assigned the same temperature 

i j Wing Angles Card 

j • SWEEPA, DIHEDA, CODEX, HMODI, TURBL (5E10.0) 

SWEEPA = wing sweep angle, in degrees 

DIHEDA = wing dihedral angle, in degrees 

CODEX = 0. for convective coefficient and adiabatic wall 

temperature computed by CAVE; = -1, for tabular 
input of coefficients and temperatures 
HMODI = nonzero value indicates that two tables will be 
read at the end and used to multiply the 
convective coefficient 

TURBL = 0. for laminar flow, = 1. for turbulent flow 
Air Properties Card (omit when CODEX = -1.) 

• GAM, RGAS, PR 

GAM = ratio of specific heats of air 
RGAS = gas constant for air, ft-lbf/lbm °R 
PR = Prandtl number of air 

Flight Trajectory Cards (omit when CODEX = -1.) 

• NTRAJ 

• TIMTAB(l), TIMTAB(2), TIMTAB(NTRAJ) 

• ALTTAB(l), ALTTAB(2), . . . , ALTTAB(NTRAJ) 

• VELTAB(l), VELTAB(2) VELTAB(NTRAJ) 

• ALPTAB(l), ALPTAB(2) ALPTAB(NTRAJ) 

NTRAJ = number of points in trajectory table (2 < NTRAJ < 50) 
TIMTAB(I) = time in trajectory table 1=1, 2, . . . , NTRAJ, sec 
ALTTAB(I) = altitude corresponding to time TIMTAB(I), ft 
VELTAB(I) = velocity corresponding to time TIMTAB(I), ft/sec 
ALPTAB(I) = angle of attack corresponding to time TIMTAB(I) must be 
non-negative, degrees 


(3E10.0) 


( 110 ) 

(8E10.0) 

(8E10.0) 

(8E10.0) 

(8E10.0) 
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Convective Coefficient and Adiabatic Wall Temperature Cards (omit when CODEX = 0.) 


The input format for these four tables is described \mder subroutine NURED 
given in Appendix D, Sheet D-1. (Note that in the structure of these tables time is 
considered argument 1 and distance argument 2, ) 

• Table of convective coefficient as a function of time and distance along 
top of leading edge 

• Table of adiabatic wall temperature as a function of time and distance 
along top of leading edge (note that the tabular entries of this table for 
time and distance must be identical to those of the above table) 

• Table of convective coefficient as a function of time and distance along 
bottom of leading edge 

• Table of adiabatic wall temperature as a function of time and distance 
along bottom of leading edge (note that the tabular entries of this table 
for time and distance must be identical to those of the above table) 

• Blank card (terminates table read in) 

Time Intervals Cards 

• NTIMES (110) 

• TIMES(l), TIMES(2) TIMES (NTIMES) (8E10.0) 

NTIMES = number of points in time intervals array (2< NTIMES <50) 
TIMES(l) = initial time (usually equals 0.), sec 
TIMES(I) = time at which temperatures will be calculated and 
printed out 1 = 2, 3, . . . , NTIMES, sec 

Convective Coefficient Modification Tables (omit when HMODI = 0.) 

Two tables are required to modify the convective heat transfer coefficient 
(see p. 14), The first table gives the multiplicative factors for the top surface of the 
leading edge; the second table gives the factors for the bottom surface. Time is 
considered argument 1 and distance argument 2. The writeup for subroutine 
NURED, Appendix D (Sheet D-1), gives the specifics on the format requirements. 
Follow these two tables with a blank card; omit the blank card if the tables are not 
read in. 
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3.3 SAMPLE PROBLEM FOR LEADING EDGE GEOMETRY 

This subsection contains an illustration of the leading edge geometry consideie 
as a sample problem (see Fig, 4). The main features are that it is made of berylliu 
has a nose radius of 0.52 cm and is cooled internally via nose and aft coolant passage 
The trajectory is one of a missile with a m^mum Mach number of six. 

Reference to the listing of the input data (see Sheets 3.3 on p. 24) shows the 
time step inteiv'als used in this problem were; 0 to 10, 10 to 20, 30 to 40.4, 40.4 to 
80, and 80 to 120 seconds. Smaller time steps were used in the beginning because 
the trajectory is changing more rapidly then. It should be noted that for a problem 
such as this one where radiation is neglected and the thermal properties are conside 
ed constant, it is possible to use one single time step to cover any portion of the 
trajectory with constant flight parameters. Specifically, since the flight parameters 
velocity, altitude, and angle of attack are constant from 40.4 to 120 seconds, it 
would have been possible to use the single interval 40.4 to 120 in lieu of the two 
intervals 40.4 to 80 and 80 to 120 seconds. For laminar flow where the convective 
coefficient is independent of wall temperature, CAVE would calculate the same 
temperatures at time 120 seconds either way since the boundary conditions and 
properties are constant throughout the interval. The 80 second point was introduced 
here so as to obtain a printout of the temperatures at this time for plotting purposes 

The following pages show listings of the input data and the output generated 
by CAVE for this leading edge problem. The sequence of the output is: 

• Statement regarding storage allocated for S array in main program 

• Geometry related input data 

• Node numbers adjacent to exterior boundary, nose cooling passage and 
aft cooling passage 

• Material properties 

• Trajectory table 

• Node number location within output array 

• Material number assigned to each node. (In this problem there is only onj 
material being used.) 



/ 

/ 


• The capacitance of each node 

• The conductance in the x-direction between nodes 

• The conductance in the y-direction between nodes 

• Initial temperature distribution 

And then the following information is printed for each time interval: 

• Flight trajectory parameters, Mach number, altitude, velocity, and 
angle of attack at the end of the time interval 

• Average heat transfer coefficients calculated using the h values at the 
beginning and end of this time interval 

• Average heat transfer couplings, which include radiation effects, if any, 
calculated using the temperatures at the beginning of this time interval 
(See Appendix C) 

• Average adiabatic wall temperatures for this time interval 

• Temperatures at the end of the interval 

• Steady-state temperatures for the heat transfer couplings and adiabatic 
wall temperatures of this interval 

• Integrated heat input to each node. This gives the net heat transfer 
at each boundary node up to the end of this time interval 

(Annotation has been added to the input and output to aid the reader.) 
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FIG. 4 GRID NETWORK* FOR MISSILE LEADING EDGE WITH COOLING 
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SHEET 3.3 OUTPUT DATA FOR LEADING EDGE PROBLEM (REFER TO FIG. 4)(SHEET 3 OF 8) 




m 

i 



.a 


29 


SHEET 3.3 OUTPUT DATA FOR LEADING EDGE PROBLEM (REFER TO FIG. 4) (SHEET 4 OF 8) 
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SHEET 3.3 OUTPUT DATA FOR LEADING EDGE PROBLEM (REFER TO FIG. 4) (SHEET 6 OF 8) 
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Section 4 


COOLED PANEL, SLAB AND X-24C STRUCTURE 
4.1 DISCUSSION 

This section discusses the input data format and definition of variables for the 
three geometries of cooled panel, slab and a structural arrangement referred to as 
basic X-24C structure. These geometries are shown in Figures 5 through 8 which give 
the overall dimensions that are input to CAVE and the type of boundary conditions 
associated with the geometries. For the cooled panel, there are two types of corners 
that can be analyzed: square and round (see Figs. 5 and 6). 

The top surface of each of the three geometries experiences aerodynamic 
heating. The user may elect to have CAVE calculate the convective heat transfer 
coefficients and adiabatic wall temperatures, or he may supply tabular inputs for 
them . If the user elects the former option, then he supplies tabular values for the 
flight parameters velocity, altitude and angle of attack as functions of time. 

In calculating the aerodynamic heating, CAVE assumes the boimdary layer 
flow is processed through an oblique shock whenever the top surface is moving at 
speeds in excess of Mach 1 and at an angle of attack relative to the freestream 
conditions. Subroutine TRANS establishes whether the flow field is laminar or 
turbulent based on a transition criterion (see Fig. B-1, Appendix B). The details 
of the aerodynamic heating equations that are used are given in Appendix B . 

Depending on the orientation of the geometry with respect to the boundary 
layer flow, the convective coefficient may vary with X (flow is from left to right in the 
plane of the paper) or be independent of X (flow is into plane of paper). Both condi- 
tions can be handled by CAVE, Setting the input variable CODEX equal to 1, selects 
the former situation, while CODEX equal to 0. selects the latter. When the convective 
coefficient is to be considered varying with X, the user must input a nonzero value for 
REFX which represents an effective boundary layer length to the left edge of the 
geometry (Refer to Appendix B) . 
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FIG. 5 COOLED PANEL GEOMETRY (SQUARE CORNER) 





AERODYNAMIC HEATING 



FIG. 7 SLAB GEOMETRY 
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As cited previously, an option exists to input directly the convective coefficient 
and the adiabatic wall temperatures as functions of time and X in lieu of the flight 
trajectory parameters. This optional feature is valuable for predicting temperatures 
within structural members subjected to other than the usual aerodynamic heating. 

The forcing functions for such special heating effects as caused by body shocks 
intersecting the wing, separated flow regions, wing-body interaction and engine 
exhaust plume impinging upon the structure can be calculated external to the code and 
then supplied as input data to the code. CODEX is set equal to -1. to exercise this 
option and the flight trajectory is not input. 

An additional option exists to modify the convective coefficients by multiplicative 
factors which are position and time dependent. (This would be useful in accounting for 
shock impingement heating.) The factors are entered as a table. A nonzero value for 
the input variable HMODI flags CAVE that this heating multiplier option will be exer- 
cised. For the normal run when the convective coefficient is not to be modified, 

HMODI equals 0. and the tables are omitted. 

This modification option can be exercised whether CAVE calculates the original 
or unmodified convective coefficients, or if they are input to CAVE. 

Subroutine SLAB2 discretizes both the cooled panel and slab geometries into 
nodes and calculates the associated capacitances and conductances. Subroutine 
X-24C discretizes the X-24C geometry. 

Figure 9 shows the grid network generated within SLAB2 for the cooled panel 
geometry of Figure 7. The network is generated based on the input dimensions SI, S2, 
Wl, W2, W3, TAU, HEIGHT, and on the AX's and AY’s (which can all be different). 
TAU represents the thickness of the insulating material which may equal zero. 

Notice that in Figure 9 nodes are located at the interface between the two materials. 
SLAB2 assumes that there are equal volumes of the two materials associated with 
each interface node. This means that one half of AY^, in this case, is associated 
with the insulator and the other half with the main material. 

Figure 10 shows the grid network generated for the cooled panel with roxmd 
comer (see Fig. 6) based on the above input dimensions plus R, the inside radius of 
the comer. 
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The convective heat transfer coefficient and the temperature of the coolant are 
input values. They are constant for the trajectory and may be input as zero, In which 
case the surface of the coolant passage is taken as being adiabatic. 

Figure 11 shows the grid network generated within SLAB2 for the slab geometry 
of Figure 7. The network is generated based on the input values for the AX's and AY’s 
(which may all be different). Aerodynamic heating takes place on the top surface and 
is handled just as in the cooled panel geometry. The other three sides of the slab are 
taken to be adiabatic. Slab geometry is handled very much the same as cooled panel 
geometry. The input card for SI, S2, Wl, etc. , associated with the dimensions of 
the cooled panel is left blank for the slab; in other cases, the required input informa- 
tion is the same. 

A review of Figure 8 shows that the X-24C geometry introduces a feature not 
embodied in any of the other geometries - contact resistance between materials. The 
X-24C geometry can be viewed as having 5 components with a total of 4 interfaces be- 
tween them. The unit surface contact resistances at the 4 interfaces are inputs to 
CAVE. Figure 12 gives the grid network generated within subroutine X-24C for the 
geometry given in Figure 8. It can be seen that nodes are located at the interfaces 
between components. The precise location of a node at an interface is in the upper 
component as the full contact resistance is applied to the conduction coupling between 
the interface node and the node below it. As with all the other geometries, the user 
must input the ifl^Y's such that the interfaces are mid-way within a AY spacing. 

Another new feature to the X-24C geometry is that up to three materials can be 
involved in the structure instead of the usual two. There is no limitation on the 
arrangement of the materials among the five components. The thermal conductivity 
and specific heats of the materials can be constant or temperature dependent. 
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AERODYNAMIC HEATING 



FIG. 11 GRID NETWORK FOR SLAB 
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4. 2 INPUT DATA FORMAT FOR COOLED PANEL AND SLAB GEOMETRIES 
Indexes Card 


• 

JGEO, L, M, NE* 


(415) 


JGEO = 0 (selects slab or cooled panel geometry) 
L = number of elements through the material 



M = number of elements along top 




NE = number of dominant eigenvalues to be used in 

solution (e. g. , a typical number is 5) 


Title Card 



• 

Rxm identification, comments, etc. 


(5A10) 

Radiation Card 



• 

EPSl, TBGl 


(2F10. 5) 


EPSl = emissivity of surface 

TBGl = background radiation temperature, 

°R 


Material Properties Cards 



• 

MAT 


(15) 

• 

NMAT, RHOl, CONAVl, CPAVl 


(I10,3F10.5) 

• 

TMATl(l), TMAT1(2), TMATl(NMATl) 

omit 

(8E10. 0) 

• 

• 

CONDTl(l), CONDTl(2), ..., CONDTl(NMATl) 
CPTl(l), CPT1(2), CPTl(NMATl) 

if 

NMAT1=0 

(8E10. 0) 
(8E10. 0) 

(If MAT = 

2 include the cards;) 



• 

NMAT2, RH02, CONAV2, CPAV2 


(I10,3F10.5) 

• 

TMAT2(1), TMAT2(2) TMAT2(NMAT2) 

omit 

(8E10. 0) 

• 

. • 

CONDT2(l), CONDT2(2), CONDT2(NMAT2) 

CPT2(1), CPT2(2), CPT2(NMAT2) 

► if 

NMAT2=0 

(8E10. 0) 
(8E10.0) 


♦Current dimension limitations require that the product L x M not exceed 200 and 
that M not exceed 50. 


MAT = number of materials (1 or 2) 

NMATl = number of entries in properties table , 

(maximum of 10) NMATl = 0 for 
constant properties 


RHOl 

CONAVl 

CP AVI 

THAT la ) 

CONDTl(I) 

CPTl(I) 


= density of material 1, Ibm/cu-ft 

- average thermal conductivity of material 1 
(used when NMATl = 0), Btu/ft-sec-TR. 

= average specific heat of material 1 (used when 
NMATl = 0), Btu/lbm-®R 

= temperatures in thermal properties table for 
which CONDTl(I) and CPTl(I) are given; 
r = 1, 2, , NMATl, °R 

= thermal conductivity of material 1 at 
temperature TMATl(I), Btu/ft-sec-°R 

- specific heat of material 1 at temperature 
TMATl(I), Btu/lbm-TEl 


NMAT2, RH02, CONAV2, etc., same as NMATl, RHOl, CONAVl, etc. 
except applied to material 2 


Detail Geometry Cards 


• DELX(l), DELX(2), DELX(3), ..., DELX(M) (8F10. 5) 

• DELY(l), DELY(2), DELY(3), ..., DELY(L) (8F10.5) 

• TAU, R (2F10.5) 

• SI, S2, Wl, W2, W3, HEIGHT, TCOOL, HCOOL . (8F10.5) 


(leave this card blank when considering a slab. Fig. 11) 


DELXa) 

DELY(I) 

TAU 

R 


spatial increments in x direction 1=1, 2, . . . , M, ft 

spatial increments in y direction 1=1, 2, L, ft 

thickness of material 1, ft 

radius of inner corner (Fig, 8) (leave blank 
if it doesn’t apply), ft (=0. for MAT=1) 
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52 
Wl 
W2 
W3 

HEIGHT 

TCOOL 

HCOOL 


= (refer to Figs. 9 or 10), ft 
= (refer to Figs . 9 or 10) , ft 
= (refer to Figs. 9 or iO) , ft 
= (refer to Figs . 9 or 10) , ft 
= (refer to Figs. 9 or lO), ft 
= (refer to Figs . 9 or’ 10) , ft 
= coolant temperature, ®R 

= coolant heat transfer coefficient, Btu/ft -sec-Tl 
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Initial Temperature Cards 


KOBE, I, T(I), II, JJ 


(215, ElO. 0,215) 


11100 

KOBE 

I 

T(I) 

n and JJ 


(15) 

= 0 or blank 
“ node number 

= node initial temperature, *R 

the node number is incremented by the spacing JJ 
“ until the limit II is reached. Each node so specified 
is assigned the same temperature 


Boundary Condition Cards 

Two options exist: (1) in the first option, the user inputs the flight trajectory 
and the code calculates the convective boundary conditions along the top surface of 
the panel in accordance with the equations presented in Appendix B; and (2) in the 
second option, the user inputs directly the convective heat transfer coefficient and 
adiabatic temperature as functions of time and distance. 

OPTION 1. FLIGHT TRAJECTORY SPECIFIED 


REFX, CODEX, HMODI 


(3E10.0) 

GAM, RGAS, PR 


(3E10.0) 

NTRAJ 


(110) 

TIMTAB(l), TIMTAB(2), .. . 

. , TIMTAB (NTRAJ) 

(8E10.0) 

ALTTAB(l), ALTTAB(2), . 

. . , ALTTAB (NTRAJ) 

(8E10.0) 

VELTAB(l), VELTAB(2), . 

. . , VELTAB (NTRAJ) 

(8E10.0) 

ALPTAB(l), ALPTAB(2), . 

. . , ALPTAB (NTRAJ) 

(8E10.0) 


REFX = effective boundary layer length, e.g. , distance from 

leading edge or nose of vehicle (refer to Appendix B), ft 
COBEX = 0. for uniform convective coefficient across top surface; 

CODEX = 1. for nonuniform convective coefficient (i.e. , 
function of x); CODEX = -1. for tabular input for 
convective coefficient and adiabatic wall temperature 
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HMODI 


= nonzero value indicates that a table will be read at the 
end and used to multiply the convective coefficients 
GAM = ratio of specific heats of air 

RGAS = gas constant for air, ft“lbf/lbm~°R 
PR = Prandtl number of air 

NTRAJ = number of points in trajectory table (2<NTRAJ<50) 
TTMTAB(I) = time in trajectory table 1 = 1, NTRAJ, sec 
ALTTAB(I) = altitude corresponding to time TIMTAB(I), ft 
VELTAB(I) = velocity corresponding to time TIMTAB(I), ft/sec 
ALPTAB(I) = angle of attack correspondiig to time TIMTAB(I), 

(must be non-negative) , degrees 

OPTION 2. CONVECTIVE COE FFICIENT AND ADIABA TIC WALL 
TEMPERATURE SPECIFIED 

• REFX, CODEX, HMODI (3E10. 0) 

REFX = effective boundary layer length, e.g. , distance from 

leading edge of nose of vehicle (refer to Appendix B), ft 
CODEX = -1. (indicates to code that Option 2 is being exercised) 

HMODI = nonzero value indicates that a table will be read at the end 
and used to multiply the convective coefficient 

Convective Coefficient and Adiabatic Wall Temperature Tables 

Two tables are required. The first gives the convective coefficient as a 

function of time (argument 1) and distance (argument 2). The second table gives 

the adiabatic wall temperature as a function of time (argument 1) and distance 

(argument 2), In setting up these tables, the same values for time and distance 

must be used in both tables. The range of the distance argument must include the 

interval REFX to REFX plus SI for the cooled panel geometry and REFX to REFX 

M 

plus WIDTH, where WIDTH = 2 for the slab geometry. The tables must be 

i=l ^ 

followed by a blank card. The specifics on the format for the tables are given in the 
descriptions of subroutine NURED in Appendix D. 
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The following input data is required for both options : 


Time Intervals Cards 

• NTIMES (110) 

• TIMES(l), TIMES(2), TIMES (NTIMES) (8E10.0) 

NTIMES(I) = number of points in time intervals array (2 ^ NTIMES < 50) 
TIMES (1) = initial time (usually equals 0. ) sec 

TIMES(I) = time at which temperature is calculated and printed 
out I = 2, 3, , NTIMES, sec 

Convective Coefficient Modification Table ( omit when HMODI = 0 . ) 

One table is required to modify the convective coefficient on the aerodynamically 
heated surface. The multiplicative factor is given as a function of time (argument 1) 
and distance (argument 2), The writeup for subroutine NURED (Sheet D-1, Appendix 
D) gives the specifics on the format requirements. Follow this table with a blank ‘ 
card; omit the blank card if the table is not read in. 



4,3 INPUT DATA FORMAT FOR X-24C GEOMETRY (See Fig. 12) 

Indexes Card 

• JGEO, L, M, NE* (415) 

JGEO = 2 (selects X-24C geometry) 

L = number of elements through the material 

M = number of elements along top 

NE = number of dominant eigenvalues to be used in solution 
(e. g. , a typical number is 5) 

Title Card 

• Run identification, comments, etc. (5A10) 

Radiation Card 

• EPSl, TBGl (2F10.5) 

EPSl = emissivity of surface 

TBGl = background radiation temperature, °R 

Material Properties Cards 


• 

MAT 


(15) 

• 

NMAT, RHOl, CONAVl, CPAVl 


(I10,3F10.5) 

• 

TMATl(l), TMAT1(2), TMATl(NMATl) 

omit 

(8E10. 0) 

• 

CONDTl(l), CONDTl(2), CONDTl(NMATl) 

CPTl(l), CPT1(2) CPTl(NMATl) 

j 

^f 

NMAT1=0 

(8E10. 0) 
(8E 10.0) 

• 

(If MAT 

= 2 include the cards :) 



• 

NMAT2, RH02, CONAV2, CPAV2 


(I10.0,3F10.5) 

• 

TMAT2(1), TMAT2(2), ..., TMAT2{NMAT2) 

omit 

(8E10. 0) 

• 

CONDT2(l), CONDT2(2) CONDT2(NMAT2) 

CPT2(1), CPT2(2), ..., CPT2(NMAT2) 

> if 

NMAT2=0 

(8E10.0) 

• 

(8E10. 0) 


MAT = number of materials (1, 2 or 3) 

NMATl = number of entries in properties table (maximum of 10) 

NMATl = 0 for constant properties 

RHOl = density of material 1, Ibm/cu-ft 

*Current dimension statements require L x M < 200 and M < 50, 
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CPAVl 


= average specific heat of material 1 (used when 
NMATl = 0), Btu/lbm-^R 
TMATl(I) = temperatures in thermal properties table for which 
CONDTl(I) and CPTl(I) are given; 1 = 1, NMATl, “R 
CONDTl(I) = thermal conductivity of material 1 at temperature 
TMATl(I), Btu/ft -sec-“R 

CPTl(I) = specific heat of material 1 at temperature 
TMATl(I), Btu/lbm-*R 

NMAT2, RH02, CONAV2, etc., same as NMATl, RHOl, CONAVl, etc., except 


applied to material 2 
Detail Geometry Cards 

• DELX(l), DELX(2), DELX(3), DELX(M) (8E10.0) 

• DELY(l), DELY(2), DELY(3), . . . DELY(L) (8E10.0) 

• SI, S2, S3, S4 (4E10.0) 

• Wl, W2, W3, W4, W5, W6, W7 (7E10.0) 

• RC(1), RC(2), RC(3), RC(4) (4E10.0) 

• MP(1), MP(2), MP(3), MP(4), M(5) (5E10.0) 


DELX(I) = spatial increments in x direction I = 1,2, . . . ,M, ft 

DELY(I) = spatial increments in y direction I = 1,2,. , . ,L, ft 

SI, S2, etc. = (Refer to Fig. 12) ft 
RC(I) = unit surface contact resistance between 

components I and I+l, sec-sq-ft-Tl/Btu 
MP(I) = integer value 1, 2 or 3 to select material properties 

for component I 

Initial Temperature Cards 


• 

• 

KODE, I, T(I), : 

# 

• 

• 

11100 
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KOBE = 0 or blank 


I = node number 

T(I) — node initial temperature , “ R 

the node number is incremented by the spacing JJ until 
n and JJ = the limited II is reached. Each node so specified is 
assigned the same temperature 


Boundary Condition Cards 


Two options exist: (1) in the first option, the user inputs the flight trajectory and 
the code calculates the convective boundary conditions along the top surface of the 
panel in accordance with the equations presented in Appendix B; and (2) in the second 
option, the user inputs directly the convective heat transfer coefficient and adiabatic 
temperature as functions of time and distance. 


OPTION 1. FLIGHT TRAJECTORY SPECIFIED 


REFX, CODEX, HMODI 


(3E10.0) 

GAM, RGAS, PR 


(3E10.0) 

NTRAJ 


(110) 

TIMTAB(l), TIMTAB(2), .... 

TIMTAB (NTRAJ) 

(8E10.0) 

ALTTAB(l), ALTTAB(2), ... 

, ALTTAB (NTRAJ) 

(8E10.0) 

VELTAB(l), VELTAB(2), ... 

, VELTAB (NTRAJ) 

(8E10.0) 

ALPTAB(l), ALPTAB(2), ... 

, ALPTAB (NTRAJ) 

(8E10.0) 


REFX = effective boundary layer length, e.g. , distance from 
leading edge or nose of vehicle (refer to 
Appendix B), ft 

CODEX = 0. for uniform convective coefficient across top 
surface; CODEX = 1. for nonuniform convective 
coefficient (i.e. , function of x); CODEX = -1, for 
tabular input for convective coefficient and adiabatic 
wall temperature 

HMODI = nonzero value indicates that a table will be read at 

the end and used to multiply the convective coefficients 

GAM = ratio of specific heats of air 
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RGAS 

PR 

NTRAJ 

TIMTAB(I) 

ALTTAB(I) 

VELTAB(I) 

ALPTAB(I) 


= gas constant for air, ft-lbf/lbm-*R 
= Prandtl number of air 

= number of points in trajectory table (2 < NTRAJ < 50) 
= time in trajectory table 1 = 1, NTRAJ, sec 
= altitude corresponding to time TIMTAB(I) , ft 
= velocity corresponding to time TIMTAB(I) , ft/sec 
= angle of attack corresponding to time TIMTAB(I) , 


in degrees 


OPTION 2. CONVECTIVE C OEFFICIENT AND ADIABATIC WALL 
TEMPERATURE SPECIFIED 

• REFX, CODEX, HMODI (3E10.0) 

REFX = effective boundary layer length, e.g. , distance from 

leading edge or nose of vehicle (refer to Appendix B) , ft 
CODEX = -1. (indicates to code that Option 2 is being exercised) 

HMODI = nonzero value indicates that a table will be read at the 


end and used to multiply the convective coefficients 


Convective Coefficient and Adiabatic Wall Temperature Tables 


Two tables are required. The first gives the convective coefficient as a ■ 

function of time (argument 1) and distance (argument 2) . The second table gives 

the adiabatic wall temperature as a function of time (argument 1) and distance 

(argument 2). In setting up these tables, the same values for time and distance 

must be used in both tables. The range of the distance argument must include the 

interval REFX to REFX plus SI for the cooled panel geometry and REFX to REFX 

M 

plus WIDTH, where WIDTH = 2 AX. for the slab geometry. The tables must be 

i=l ^ 

followed by a blank card. The specifics on the format for the tables are given in 
the descriptions of subroutine NURED in Appendix D. 

The following input data is required for both options: 

Time Intervals Cards 

• N TIMES (110) 

• TIMES(l), TIMES(2) TIMES (NTIMES) (8E10.0) 

NTIMES(I) = number of points in time intervals array (2 < NTIMES < 50) 
TIMES(l) = initial time (usually equals 0.) sec 

TIMES(I) = time at which temperature is calculated and printed 
out 1 = 2,3,... , NTIMES, sec 
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Convective Coefficient Modification Table (omit when HMODI = 0.) 


One table is required to modify the convective coefficient on the aerodynamically 
heated surface. The multiplicative factor is given as a function of time (argument 1) 
and distance (argument 2). The writeup for subroutine NURED (Sheet D-1, Appendix 
P) gives the specifics on the format requirements. Follow this table with a blank 
card; omit the blank card if the table isn’t read in. 
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4.4 SAMPLE PROBLEM FOR COOLED PANEL GEOMETRY 

This subsection contains an illustration of the cooled panel geometry considered 
as a sample problem (see Sheet 4.4). Figures 9 and 10 show the grid network. The 
panel was constructed of aluminum with a 0.686 cm layer of beryllium applied across 
the tpp in contact with the hot boundary layer air. The aluminum structure was cooled 
by a coolant at a temperature of 660°R with a convective heat transfer coefficient of 
612.8 watts/m^K. 

Aerodynamic heating was calculated by CAVE based on the missile flight trajec- 
tory shown in the following listing of the input data. The panel was located on the mis- 
sile 6.096 m aft of the start of the boundary layer and hence the input variable 
REFX equals 20. Throughout the flight trajectory the angle of attack was 20 degrees. 

This sample problem shows the use of the heating multiplier option. Examina- 
tion of the heating multiplier table in the listing of the input data reveals that during 
the time period 30 to 40.4 seconds, the convective heat transfer coefficient increases 
fivefold for the center section of the panel. To represent this step change in the 
heating, table entries of 29.9, 30. 1, 40.3 and 40.5 have been used. Furthermore 
the time step intervals array contains these same times. Whenever a step change 
is to be approximated to CAVE, it is necessary to use ’’bracketing” times in setting 
up a multiplier table and the time step interval array. That is, if a step change 
occurs at time t^, the table should have time entries at t^- € and t^+ € where € is 
some small number. Likewise the time step interval array should have these same 
times. 

The following pages show listings of the input data and the output generated by 
CAVE for this cooled panel problem. The sequence of the output is:. 

• Statement regarding storage allocated for S array in main program 

• Geometry related input data 

• Material properties illustrating the format used when the properties are 
temperature dependent and temperature independent 

• Flight trajectory table 

• Table of heating multiplier as a function of time and distance 
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• Node number location within the output array . A node number of 0 (zero) 
printed in this array indicates a nonactive node . Notice how the printout has 

a format similar to the shape of the geometry of the problem being considered, 
which makes for very convenient reading and quickly gives the user a feel 
for the temperature distribution within the body. This format is carried 
throughout in all the following arrays 

• Material number assigned to each node, A material number of 3 signifies 
that the node is at an interface between materi al 1 and 2 

• Capacitance of each node 

• Conductance in the x-direction between nodes 

• Conductance in the y-direction between nodes 

• Initial temperature distribution 

And then , the following information is printed for each time interval ; 

• Flight trajectory parameters, Mach number, altitude, velocity, and angle of 
attack at the end of the interval 

• Average heat transfer coefficients calculated using h values at the beginning 
and end of this time interval 

• Average heat transfer couplings, which include radiation effects, if any, 
calculated using the temperatures at the beginning of this time interval 
(See Appendix C) 

• Average adiabatic wall temperatures for this time interval 

• Temperatures at the end of the interval 

• Steady-state temperatures for the heat transfer couplings and adiabatic wall 
temperatures of this interval 

• Integrated heat input to each node. This gives the net heat transfer at 
each boundary node up to the end of this time interval, A positive value 
signifies heat transfer into the body and a negative value heat transfer out 
of the body at that boundary node 

(To aid the reader, annotation has been added to the listing of the input data and 
to the output . ) 
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SHEET 4.4 INPUT DATA FOR COOLED PANEL PROBLEM PREFER TO FIGS. 9 & 10) (SHEET 1 OF 8) 
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SHEET 4.4 OUTPUT DATA FOR COOLED PANEL PROBLEM (REFER TO FIGS. 9 & 10) (SHEET 2 OF 8) 
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SHEET 4.4 OUTPUT DATA FOR COOLED PANEL PROBLEM (REFER TO FIGS. 9 & 10) (SHEET 3 OF 8) 
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SHEET 4.4 OUTPUT DATA FOR COOLED PANEL PROBLEM (REFER TO FIGS. 9 & 10) (SHEET 4 OF 8) 






cnnmJCT»Ncc in v-niREcxioN 
iw / cdl I 



SHEET 4.4 OUTPUT DATA FOR COOLED PANEL PROBLEM (REFER TO FIGS. 9 & 10) (SHEET 5 OF 8) 
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SHEET 4.4 OUTPUT DATA FOR COOLED PANEL PROBLEM (REFER TO FIGS. 9 & 10) (SHEET 7 OF 8) 








SHEET 4.4 OUTPUT DATA FOR COOLED PANEL PROBLEM (REFER TO FIGS. 9 & 10)(SHEET 8 OF 8) 







4. 5 SAMPLE PROBLEM FOR SLAB GEOMETRY 


Figure 11 illustrates the grid network for the slab geometry considered as a 
sample problem. Aluminum was used throughout the slab and hence TAU = 0. 

This sample problem illustrates the use of the option to input the convective 
heat transfer coefficients, h, and adiabatic wall temperatures, as functions of 

time and distance along the top surface rather than have CAVE calculate them based 
on a flight trajectory. The values for h that were input were based on the equation 

h = 0.005 + 0.0075 cos(5 7Tx) 

The adiabatic wall temperature was taken to be constant with respect to both time 
and distance. 

Since the slab is rectangular in shape , the output is particularly easy to read; 
a row in an output array giving the values associated with a row of nodes. The type 
of output is the same as was described in detail for the cooled panel geometry 
(see subsection 4.4). 
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SHEET 4.5 INPUT DATA FOR SLAB PROBLEM (REFER TO FIG. 11)(SHEET 1 OF 6) 
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SHEET 4.5 OUTPUT DATA FOR SLAB PROBLEM (REFER TO FIG. 11) (SHEET 2 OF 6) 
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SHEET 4.5 OUTPUT DATA FOR SLAB PROBLEM (REFER TO FIG. 11) (SHEET 3 OF 6) 
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4. 6 SAMPLE PROBLEM X-24C GEOMETRY 

Figure 12 shows the grid network for the X-24C geometry considered as a 
sample problem . 

In many respects the input and output for this problem are similar to that of the 
cooled panel problem presented in subsection 4.4. One of the differences is that 
there are three materials here instead of two which affects the array of material 
.numbers. 

For the X-24C geometry a material number of 1, 2 or 3 at a node signifies that a 
node is comprised entirely of material 1, 2 or 3. (Recall that the properties of the 
materials is established via the input. ) And a two-digit material number signifies that 
the node is at an interface between two materials, with the tens digit giving the 
material number for the upper material and the ones digit the material number for 
the lower material at the interface. For example, node 54 is located at the inter- 
face between the beryllium and the insulation, since beryllium has been set up to be 
material 1 and the insulation to be material 2 we find the material number of node 54 
to be 12. Material numbers are set up by subroutine X-24C and they provide a check 
on whether geometry data has been input correctly to CAVE. (The input and output 
listings on the following pages have been annotated. ) 
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Section 5 


GENERAL GEOMETRY 


5. 1 DISCUSSION 


This section describes the general geometry capabilities of CAVE and gives the 
input data format. 


When exercising the general geometry option of CAVE , the user must do the 
usual tedious and laborious calculations associated with setting up a thermal network. 
Namely, the volume of each node, the cross-sectional area divided by the X distance 
between adjacent nodes, and the cross-sectional area divided by the Y distance 
between adjacent nodes all must be supplied as input data to the code. The code will 
then multiply the volumes by p Cp to obtain the node capacitances . The area over 
distances will be multiplied by the thermal conductivity to obtain the conductances . 
Multimaterial problems can easily be handled by supplying capacitances instead of 
volumes and conductances instead of area over distances as input data. The material 
properties p , and k should be input as 1 in this case. 


It is possible to simulate within CAVE a convection coupling, either constant 
or time varying , to each node . The values of the couplings are supplied as input 
data. Radiation heat transfer is not considered by CAVE when the general geometry 
option is selected. 


The matrix package within CAVE is very efficient from both an execution and 
storage standpoint. This is made possible by some limitations that have been incor- 
porated into the package. Referring to Figure 13, the limitations may be stated 
as follows : 


1 . 


There exists a conductance coupling between node I and node I + 1, e.g. , 
nodes 16 and 17. (The value of the coupling may be zero.) One exception 
is when node I is in the bottom row. For example, nodes 12 and 13 are 
not coupled . 
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L = 6, M = 7, N = 6x7 = 42 


FIG. 13 GRID NETWORK FOR A GENERAL GEOMETRY PROBLEM 
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2 . There exists a conductance coupling between node I airi node I + L where 
L is the number of elements in the Y direction, e.g. , nodes 16 and 22. 

(The value of the coupling may be zero . ) One exception to this is when 
node I is in the right column. Since there is no node I + L, there can 
be no conductance coupling. 

3 . Node numbering is done sequentially starting with 1 and going to n , the 
number of nodes . Numbers run columnwise starting at the top of the left 
column with 1, the node below it 2 and the node to the right of it 1 + L. 

4 . No other conductance couplings exist . That is , node I will have at most 
couplings with nodes I-l, I+l, I-L and I+L. It cannot be coupled to any 
other nodes . 

The last assumption precludes, for example, a nodal network having the couplings 
shown in Figure 14. 

On occasion, fictitious nodes are present within the network. They are forced 
into the network by virtue of the above assumptions regarding the couplings. No input 
data is required for these nodes since they are not ar active part of the problem . 

Nodes 3, 4, 5, 19, 20, 21 and 27 (given in Fig. 13) are such fictitious nodes. 

To visualize the class of problems that CAVE can most readily accommodate, 
consider a rectangular grid network that gets cut up to form different shapes . If the 
grid network is drawn on a rubber sheet and that sheet is stretched at will , any 
resulting figure could be analyzed by CAVE. 

The following subsections give the input data format and definition of input 
variables for running a general problem . 
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FIG. 14 PORTION OF A GRID NETWORK THAT CANNOT BE HANDLED 


I 

f ' 

I ' 

Ij 5.2 INPUT DATA FORMAT FOR GENERAL GEOMETRY (Fig. 13) 

Basic Geometry Card 

’ • JGEO, L, M, NE* (415) 

L 

I JGEO = -1 (selects general geometry option) 

L = maximum number of elements in the 
Y direction 

M = maximum number of elements in the 
X direction 

NE - number of dominant eigenvalues to be used 
in solution (e.g. » a typical number is 3) 

Title Card 

• Run identification, comments, etc. (5A10) 

• Blank Card 

Material Properties Cards 


. • 

MAT 



(L5) 

• 

NMATl, RHOl, CONAVl, 

CPAVl 


(110, 3F10.5) 

• 

TMATl(l), TMAT1(2), ... 

. , TMATl(NMATl) 

omit 

(8E10.0) 

• 

CONDTl(l), CONDTl(2), 

. . . , CONDTl(NMATl) 

^if 

(8E10.0) 

• 

CPTl(l), CPT1(2), ..., CPTl(NMATl) 

NMAT1=0 

(8E10.0) 

(If MAT 

= 2 include the cards ;) 




• 

NMAT2, RH02, CONAV2 , 

CPAV2 


(I10.0,3F10.5) 

• 

TMAT2 (1) , TMAT2 (2) , . . 

. , TMAT2(NMAT2) 

omit 

(8E10.0) 

• 

CONDT2(l), CONDT2(2), 

. . . , CONDT2(NMAT2) 

^if 

(8E10.0) 

• 

CPT2(1), CPT2(2), ..., CPT2(NMAT2) I 

NMAT=2 

(8E10.0) 


MAT = number of materials (1, 2 or 3) 

NMATl = number of entries in properties table (maximum of 10) 

NMATl = 0 for constant properties 
RHOl = density of material 1, Ibm/cu-ft 

CONAVl = average thermal conductivity of material 1 (used when 
NMATl = 0) , Btu/ft-sec-“R 

CPAVl = average specific heat of material 1 (used when 
NMATl = 0), Btu/lbm-“R 

^Current dimension limitations require that the product LXM not exceed 200. 
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/ 

I 

TMATl(I) = temperatures in thermal properties table for which 

CONDTl<I) and CPTl(I) are given; 1=1, 2, . . . , NMATl, °R 
CONDTl(I) = thermal conductivity of material 1 at temperature 
TMATl(I), Btu/ft-sec-®R 

CPTl(I) = specific heat of material 1 at temperature THAT 1(1) , 

Btu/lbm-“R 

NMAT2, RH02, CONAV2, etc. , same as NMATl, 

RHOl, CONAVl, etc. , except applied to material 2 

Volume Cards 

• KOBE, I, V(I), II, JJ (215, E 10. 0,215) 

w • • • 

V • • • 

V • • • 

• 11100 (15) 

KOBE = 0 or blank 

I = node number 

V(I) = node volume, cu ft (or optionally pVC ), Btu/°R 

P 

n = limit for multiple parameter input 
JJ = increment for multiple parameter input 

(The node number is incremented by the spacing JJ until 
the limit H is reached. Each node so specified is assigned 
the same temperature . ) 

Material Selection Cards 

• KOBE, I, MATNUM(I), R, JJ (215 ,E10. 0,215) 

• 11100 
KOBE 
I 

MATNUM(I) 


as) 

= 0 or blank 
= node number 

= 1 to use material 1 properties 
= 2 to use material 2 properties 

= 3 when node is at interface between materials 1 and 2 
(interface is understood to be parallel to X-axis) 
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II 

JJ 


= limit for multiple parameter input 
= increment for multiple paremeter input 


Area Over X Cards 

• KODE, I, AOVERX(I), n, JJ 


(215, E 10. 0,215) 


moo 

KODE 

I 

AOVERX(I) 


n 

JJ 


(15) 

0 or blank 
node number 

cross-sectional area divided by X distance between 
node I and node I + L, required only for nonzero 
conductances, ft2/ft (or optionally, the conductance 
between nodes I and I + L) , Btu/ sec-“R 
limit for multiple parameter input 
increment for multiple parameter input 


Area Over Y Cards 


KODE, I, AOVERY(I), II, JJ 


(215, ElO. 0,215) 


11100 

KODE 

I 

AOVERX(I) 


n 

JJ 


(15) 

0 or blank 
node number 

cross-sectional area divided by Y distance between 
node I and node I + 1, required only for nonzero 
conductances, ft^/ft (or optionally, the conductance 
between nodes I and I + 1) , Btu/sec-°R 
limit for multiple parameter input 
increment for multiple parameter input 
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Convection Coupling Cards 


(215, ElO. 0,215) / 


KOBE, I, HA(I), II, JJ 


11100 

KOBE 

I 

HA(I) 


II 

JJ 


(15) 

= 0 or blank 
= node number 

= convective coupling between node I and TAW(I) , 
Btu/sec-“R 

(or optionally, a negative number to indicate that 
the coupling is time varying. The absolute value 
of the negative number is then the number of the 
dependent variable in the convection table that 
gives HA (I) as a function of time) 

= limit for multiple parameter input 
= increment for multiple parameter input 


Adiabatic WaU Temperature Cards 


KOBE, I, TAW (I), II, JJ 


(215, ElO. 0,215) 


11100 

KOBE 

I 

TAW(I) 


n 

JJ 


(15) 

= 0 or blank 
= node number 

= adiabatic wall temperature associated with node I, °R 
(or optionally, a negative number to indicate that the 
adiabatic wall temperature is time varying. The absolute 
value of the negative number is then the number of the 
dependent variable in the convection table that gives TAW 
as a function of time) 

= limit for multiple parameter input 
= increment for multiple parameter input 
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Convection Table Cards 


• 

KODE, LI, L2, TITLE 

(315, A65) 

• 

TIME (I), TIME (2), . . . , TIME (LI) 

(7E10.0) 

• 

Dl(l), Dl(2), ... , Dl(Ll) 

(7E10.0) 

• 

• 

A 

D2(l), D2(2) D2(L1) 

(7E10.0) 

w 

• 

• 

DL2(1), DL2(2), ..., DL2(L1) 

(7E10.0) 

• 

11100 

KODE = 0 or blank 

LI = number of values of time , the independent 

variable (LI < 250/L2) 

L2 = number of dependent variables (1 < L2 < 20) 

(15) 


TIME (I) = time in convection table 1=1,2, ... LI, sec 

D1(I) = first dependent variable to represent either 

HA or TAW associated with a node or nodes 
at TIME (I) 

• ... 

A 


V • • • 

• ... 

DL2 (I) = last dependent variable to represent either HA or TAW 

associated with a node or nodes at TIME (I) 

Initial Temperature Cards 


• 

• 

A 

KODE, I, T(I), II, JJ 

(215, E 10. 0,215) 

• 

• 

11100 

(15) 


KODE = 0 or blank 

I - node number 
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T(I) 

n and JJ 


Time Interval Cards 


/ 

= node initial temperature , ®R 

! the node number is incremented by the spacing JJ 
imtil the limit II is reached . Each node so specified 
is assigned the same temperature 




NTIMES (110) 

TIMES(I), TIMES(2) TIMES(NTIMES) (8E10.0) 


NTIMES = number of points in time intervals 
array (2 < NTIMES < 50) 

TIMES (I) = initial time (usually equals 0.)f sec 

TIMES (I) = time at which temperatures will be calculated 

and printed out 1 = 2, 3 , . . . , NTIMES , sec 


5.3 SAMPLE PROBLEM FOR GENERAL GEOMETRY (Fig. 15) 


The general capabilities of CAVE are illustrated via the simple grid network 
shown in Figure 15. While CAVE can handle much more arbitarily shaped geometries 
than this one, all the basic ingredients are present here. 

As Figure 15 shows , two materials were used: a layer of beryllium across the 
top that was being aerodynamic ally heated and the remainder aluminum with a portion 
of its boundary being cooled. The convective coefficients and adiabatic wall temper- 
atures associated with the aerodynamic heating are shown as functions of distance and 
time in Figure 15. The values of hAX and for the different times at nodes 1, 9, 
17, 25 and 33 were provided as input data to represent the aerodynamic heating. 

The following pages present listings of the input data used for this problem and 
the resulting output data. The input data for this general type of geometry is seen 
to be more extensive than for one of the built-in configurations. This is due to the 
need to input for each node a volume, material number, condition area over distanceSt 
convection coupling and adiabatic wall temperature, in addition to the usual material 
properties initial temperatures and time step intervals. 

The output from CAVE for a general geometry problem is very similar to that 
of a built-in geometry. At the end of every time interval there is a printout of the 
convection couplings, adiabatic wall temperatures, steady-state temperatures and 
integrated heat inputs. No printout of convective coefficients appears since the 
general geometry option deals strictly with the coupling, hAx and not with h alone. 

(The input and output have been annotated to assist the reader.) 
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SHEET 5.3 INPUT DATA FOR GENERAL GEOMETRY PROBLEM (REFER TO FIG. 15) (SHEET 2 OF 6) 
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SHEET 5.3 OUTPUT DATA FOR GENERAL GEOMETRY PROBLEM (REFER TO FIG. 15) (SHEET 3 OF 6) 
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APPENDIX A 


DESCRIPTION OF THE HYBRID ANALYTICAL-NUMERICAL TECHNIQUE. 

This appendix presents a brief summary of the Hybrid Analytical-Numerical 
(HAN) technique using, as an example, a one -dimensional, finite thickness, conduction 
problem. This application provides two services. First it provides a clear exposition 
of the HAN technique with its attendant matrix operations ; and second , it provides 
insight into the accuracy of the. technique since an exact solution to the problem is 
known. The effect of retaining only the dominant eigenvalues and eigenvectors (E & E) 
on the accuracy of the solution can be clearly assessed. In NASA CR-2435 , Maise 
and Rossi thoroughly investigated the effect that the number of E & E’s have on the 
accuracy of the h predicted for the inverse problem. In addition, for a direct prob- 
lem at a particular time , they show the typical ei jors incurred by neglecting 
subdominant eigenvalues in a very simple model problem. This appendix 
explores how the incurred errors vary as a function of time . The solution to the 
sample problem was obtained using a specially prepared computer code to perform 
all of the matrix operations and companion calculations . The solution provided a 
completely independent check on the eigenvalues , eigenvectors and temperatures 
generated by CAVE. 

The problem considered was that of a slab heated by convection on one face and 
perfectly insulated on the other face. Figure A-1 shows the 10-node network selected 
to solve the problem; also shown are the properties, temperatures and dimensions 
used. 

In the HAN method , spatial derivatives are replaced by their appropriate finite 
difference representations and the temporal derivatives are retained as ordinary 
derivatives. Ri effect, the problem is subdivided into a number of uniform tempera- 
ture systems or nodes that are coupled and changing in temperature. Utilizing the 
notation of Appendix F, the set of ten, first-order, linear, ordinary differential 
equations for the temperatures at the ten nodes can be written as; 

MT = BT + F Eq. (A-la) 
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subject to the initial condition, 

r 100 '1 


T (0) = = 





Eq. (A -lb) 


L 100 J 


The matrix is a 10 x 10 diagonal matrix associated with the heat capacity 
of the nodes and is given by; 


1.65 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

3.30 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0,0 

0.0 

0.0 

3.30 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

3.30 

0.0 

0.0 

0.0 

0.0 

0.0 

0,0 

0.0 

0.0 

0.0 

0.0 

3 . 30 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

3.30 
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0.0 
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0.0 

0.0 

0.0 
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3. 30 
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0.0 

0.0 
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0.0 

0.0 

0.0 

0.0 

3 . 30 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

1,65 


where , 


= M 


Ax 


10 10 


P c 


M.. = AxPC for i = 2, 3, 


The matrix B, also 10 x 10 in size, is associated with the heat transfer 
couplings between nodes and is given by: 



- 2000 . 
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0 . 
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0 . 

0 . 
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where , 


'11 


= -h - 


Ax 


b.. 

II 


= -2 


Ax 


= 

*10 10 Ax 


and, 


b. . , , = b. 


i,i+l 1+1, i Ax 


i =2, 3, . 


9 


i= 1, 2, ... 9 


The ten-component column vector F, represents the forcing function in the 
problem. It is given by: 


' 000 
0 
0 
0 


0 

0 

0 

0 


The ten-component volume vector T represents the temperature at the ten 
nodes at any instant of time. The vector T is the time derivative of T. 

From Appendix F, the analytic solution to the initial value problem (A-1) is 
given by: 

T = V exp (At) (T.^.^ - T^ ) Eq. (A-2 


F = 


2C 


where , f ^ = h T 


AW 
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where 


= the steady-state solution to (A-1) 

A = a diagonal matrix formed with the eigenvalues X . of matrix A 
V = a matrix of eigenvectors of matrix A 

and, 

A = defined by A = 

Specifically, these quantities are given by: 

r 200 ■) 



200 


- 7.48 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

- 67.34 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 - 

185.27 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 - 

352.85 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 - 

553.36 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 - 

764.06 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0,0 

0,0 - 

959 . 34 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0 . 0 - 

1113.2 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0 . 0 - 

1200.5 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0 . 0 - 

1463.2 
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0.221 

- 0.235 
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- 0.841 

- 0.138 

- 0.357 

0.455 

0.444 

- 0.361 

- 0.244 

- 0.126 

- 0.038 

0.002 

0.493 

- 0.203 

- 0.445 

0.363 

0.058 

0.269 

0.457 

0.451 

0.299 

- 0.099 

- 0.204 

- 0.263 

- 0.435 

0.049 

- 0.396 

0.408 

0,006 

- 0.399 

- 0.462 

0.192 

0.085 

- 0.316 

- 0. 328 

- 0.295 

- 0.388 

- 0,198 

- 0.460 

0.014 

0.475 

- 0.278 

- 0.035 

- 0.362 

- 0.148 

- 0.459 

0.071 

- 0.442 

0.234 

0.382 

- 0.332 

0.353 

0.015 

- 0.399 

0.065 

- 0.342 

0.448 

0.121 

0.338 

- 0.460 

0.081 

- 0.415 

- 0.006 

- 0.426 

0.263 

- 0.017 

0.303 

0.463 

- 0.410 

0,154 

0.196 

0.461 

0.003 

- 0.442 

0.403 

0.319 

- 0.194 

- 0.041 

- 0.124 

0.280 

- 0.409 

- 0.489 

- 0.001 

- 0.317 

0.321 

0.325 

- 0.329 

- 0.333 

0.336 

- 0.340 

0.346 

0.352 

0.001 


where the columns of V are eigenvectors of matrix A. 
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For coxnputationa.1 purposes the solution (A-2) is best written as: 


T = + C P (t) 


Eq. (A-3) 


where the matrix C is given by 
C = V D 


Here D is a diagonal matrix constructed from the components of the column 
vector R, defined by: 




(I 


init 


-I«) 


That is , the diagonal elements of D are given by: 


d 


ii 


r 


i 


The time dependent column vector P is given by: 


P(t) = 


exp (-7.48t) " 

exp (-67.34t) 
exp (-185.3t) 
exp (-352. 9t) 
exp (-553. 4t) 
exp (-764. It) r 
exp (-959. 3t) 
exp (-1113. t) 
exp (-1201. t) 
exp (-1463.t) 


That is , the i*^ 


component of P is given by : 


Pj = exp (Ajt) 
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The matrix C for the particular problem being considered is given by: 


- 1 * 9.47 

- 16.04 

- 11.80 

- 8.42 

- 6.03 

- 4.32 

- 2.90 

- 1.50 

- 0.23 

- 29.29 

- 38.70 

- 30.29 

- 20.00 

- 11.93 

- 6.56 

- 3.19 

- 1.21 

- 0.25 

- 0.00 

12.13 

- 56.98 

- 37.82 

- 15.96 

- 1.55 

4.89 

5.98 

4.32 

1.91 

0.24 

- 5.03 

- 73. 85 

- 36.94 

- 2.17 

10.63 

7.41 

0.07 

- 3.82 

- 2.96 

- 0.46 

2.08 

- 88.89 

- 27.85 

12.95 

10.44 

- 3.60 

- 6.02 

0.14 

3.03 

0.67 

- 0.86 

- 101.74 

- 12.58 

20.15 

- 1.91 

- 8.03 

3.07 

3.66 

- 2.12 

- 6.85 

0.36 

- 112.08 

5.50 

15.03 

- 12.04 

2.21 

4.42 

- 4.41 

0.52 

1.00 

- 0.15 

- 119.66 

22.34 

0.73 

- 8.15 

8.42 

- 5.37 

1.48 

1.25 

- 1.11 

0.06 

- 124.28 

34.23 

- 14.03 

5.23 

- 0.74 

- 1.62 

2.68 

- 2.62 

1.17 

- 0.03 

- 125.83 

38.51 

- 20.20 

12.52 

- 8.55 

6.21 

- 4.61 

3.13 

- 1.20 

0.02 


th 

Returning to Eq. (A -3) it can be seen that the temperature at the i node is 
given by: 

N 

T. =T^. + ^ C-. exp (X.t) Eq. (A-4) 

j=l 


where N equals the number of nodes (10 in this particular case) . 
For example, the temperature at node 1 is given by: 


T^ = 200 - 19.47 exp (-7.48t) - 16.04 exp <'-67.34t) - 11.80 exp (-185. 35t) - i.. 

Eq. (A-5) 


while the temperature at node 10 is : 


Tj^q = 200 - 125.8 exp (-7.48t) + 38.51 exp (-67.34t) - 20.20 exp (-185. 35t) + . .. 

Eq. (A-6) 

An examination of the arguments of the exponential functions shows that the 
leading terms in the finite series dominate the sum. 

Two questions naturally arise: (1) how does the HAN solution given in Eq. (A-4) 
compare, both in form and accuracy, with the exact analytical solution; and (2) what 
loss in accuracy is incurred by retaining only the dominant E & E^s in the solution 
(i.e. , the first few terms in Equation (A-4). These questions are examined in the 
foilo'wing discussions. 
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The analytical solution to the problem being considered is given by the 
solution to the partial differential equation: 


afr 1 ai 


subject to the boundary conditions , 




^ (L,t) = 0 


with initial condition , 

T0^.O) = Tj^lt 

The solution to this boundary value problem can be obtained by separation of 
variables in conjunction with Fourier series expansions or simply by looking up the 
solution to a similar problem in Car slaw and Jaeger ^s, "Conduction of Heat in 
Solids" on page 122, In either case, the solution is given by: 


M 2H sec (Aj) 2,2 

exp(-^^o.t/L > 


where the m j satisfy the equation 


tan /i . - H (see Table 1 of Carslaw and Jaeger, p. 491) 


and H =■ 
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Eq. (A*-7) shows that for a given value of x the solution is the form 


Ti(t)=TA^+ ^ a.. exp(b.t) Eq. (A-8) 

where the i subscript is used to connote that the temperature is for a particular value 
of X and where 


and 


^ij ■" ^init “ 


2H sec(A.) 
^ 3 



Eq. (A-8) is identical in form with Eq. (A -4) except that in Eq. (A-8) the series 
is an infinite one rather than a finite one. (Note that in Eq. (A-4) represents the 
steady-state solution at node i , which is in fact • 

It is instructive to compare the C^^^s with the a^.’s and the X.'s "with the b/s. 
These comparisons are shown in the following tables for the two ejitreme locations 
for X, X. e. , X - 0 and x = L corresponding respectively to i = 1 and i = 10. 



1 = 1 

X 

It 

o 

i = 10 

"■■■ 

x= L 

j 

<=ii 


^lOj 

®10j 

1 

-19.47 

-19.56 

-125.83 

-125.98 

2 

-16.04 

-16.63 

38.51 

38.80 

3 

-11.80 

-12.72 

-20.20 

-20.40 

4 

-8.42 

-9.33 

12.52 

12.46 

5 

-6.03 

-6.84 

-8.55 

-8.28 

6 

-4.32 

-5.10 

6.21 

5.83 
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The X /s and b.’s are- ind^ endent of x. 


They are as follows: 


i 



1 

-7.48 

-7.49 

2 

-67.34 

-68.19 

3 

-185.27 

-192.90 

4 

-352.85 

-385.43 

5 

-553.36 

-648.50 

6 ■ 

-764.06 

-986.60 


The comparison between the sets of constants is very good. This is to be 

expected since it is known that as the number of nodal points N approaches infinity the 

C..^s will approach the a./s and the X/s will approach the b.’s, and, therefore, in 
1] 13 3 3 

the limit, the HAN solution approaches the exact solution. But the really pertinent 

question here is: how accurate is the 10-point HAN solution? Figures A-2 and A-3 
address this question. They show comparisons of the HAN solution with the exact 
solvttion at the front and rear faces of the slab. If all ten terms (i.e. , eigenvalues 
and eigenvectors) of the HAN solution are used , the temperatures calculated are 
virtually identical to those of the exact solution. As a matter of fact for times in 
excess of 0. 05 hr , the first exponential term in the HAN solution compares with the 
exact solution to within 0.06® 1C It is very interesting to observe that the shorter the 
time period of interest the greater the number of terms required to achieve a given 
accuracy. An alternative point of view to this observation is perhaps more significant. 
Namely, for a given number of terms or eigenvalues there is a minimum time period 
to achieve a certain accuracy. It may be noted that (Ax) /a can be considered a 
characteristic time for the transient response of nodal elements. For this example, 
(Ax) / a cs 0.003 hr, and it can be seen that the time required for accuracy with only 
a few dominant eigenvalues is of the order of three times this value. A further note is 
that the typical explicit finite algorithm has a maximum time step (for the problem 
being studied it is 0.00165 hr). This maximum time step arises from stability require- 
ments, and is not related to the accuracy considerations for the HAN method. 

Two things need emphasizing at this juncture. First, with the HAN method con- 
siderable computer machine time can be saved by calculating only those E & E’s that 
are "significant". This was noted very aptly by Maise and Rossi in NASA CR-2435 and 
used by them in the CAPE code for the indirect heat transfer problem. As regards the 
direct heat transfer problem being solved here, Figures A-2 and A-3 show that for a 
time period of interest greater than 0.01 hr, it would be a waste of machine time to 
compute the fourth through tenth E & E’s since they would have an insignificant effect 
on the computed temperature. 
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FIG. A-3 COMPARISON OF NODAL POINT EIGENVALUE SOLUTION AND THE EXACT SOLUTION AT X = 0.9 





The second item that needs emphasizing is that because the boundary conditions 
on the problem studied herein were constant with respect to time , it was necessary to 
obtain the E & E*s only once. Time var3dng boundary conditions will require periodic 
revision of the E & E's to reflect the changes in matrix A. In essence, problems 
involving time varying boundary conditions are solved by dividing the time interval 
into M subintervals within which the boundary conditions will be constant. Thus there 
will be M subproblems with the initial conditions of one problem being the final" 
conditions of the previous problem. Obviously, frequent revisions imply a small 
time period and hence more E & E*s to achieve a given accuracy. Some reflection on 
all of this reveals that the optimal solution (viz. , minimum machine time for a given 
accuracy level or maximum accuracy level for a given machine time) will be achieved 
by proper selection of tie number of E & E revisions (or time steps) and the number 
of E & E’s to be calculated. The following table illustrates this point. In Table A-1, 
M represents the optimal number of E & E's revisions and NE represents the optimal 
number of E & E's. 


Table A-1. Mix of Constant Machine Time Solutions 


Time Period 
Between E & E 
Revision 

No. of 
E & E 
Revisions 

No. of 
E & E's 

Remarks 

Long 

<M 

>NE 

With long time period, do not need this 
many E & E's. Furthermore, the piecewise 
constant representation of boundary 
condition is poor 

Optimal 

M 

NE 

Number of E & E's calculated consonant 
with the time period. Good piecewise 
constant representation of time varying 
boundary conditions 

Short 

>M 

<NE 

Number of E & E's insufficient to provide 
good accuracy for this short time jDeriod. 
Piecewise constant representation of 
boundary condition is excellent 



The '’optimal" values for M and NE depend somewhat on the particular flight 
trajectory, geometry, dimensions and materials being analyzed. For the missile 
and Hypersonic Research Aircraft problems that were analyzed during the develop- 
ment of the CAVE code , NE should be approximately 3 to 5 , and M should be 
approximately twice the number of flight trajectory points input to the code . 

In summary, this appendix has presented the details of the HAN technique 
applied to a one -dimensional conduction problem. The HAN solution gave a near 
perfect comparison with the exact analytical solution of the same problem. It was 
seen that the HAN solution has a minimum time interval for a given accuracy level 
and number of E & E's. Furthermore, we saw that there is an optimal number of 
time subdivisions and E & E's for a specified accuracy level or machine time. 
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APPENDIX B 


AERODYNAMIC HEATING EQUATIONS* 

Incorporated within the CAVE code are equations to predict the aerodynamic 
heat transfer that will occur at the surfaces of the leading edge and flat plate geome- 
tries. This appendix presents the equations that are used in subroutines OVLY20, 
LEESl and FLATH. The equations are valid up to a Mach number of six or so, where 
real gas effects become important. Ambient pressure and temperature as functions 
of altitude are calculated in subroutine ATMOS based on the 1962 U.S, Standard Atmos- 
phere. The range of altitude is from sea level to 47 350 m. 



Aerodynamic heating of the leading -edge geometry is handled in two steps : 

(1) the convective heat transfer coefficient at the stagnation point is computed in sub- 
routine OVLY20 and; (2) the ratios of the local convective values to the stagnation 
point value are calculated in subroutine LEESl for all the surface nodes. 


The user can flag CAVE to use either a turbulent or laminar flow correlation 
for the stagnation point coefficient. 

For turbulent flow, CAVE uses Beckwith and Gallagher's equation (B-1); 


2R 


Pr (o . 022 8 


N 


T D \ 

Mo W 


|376 ^00 


cos A 


eff 


2Rn 

dVJ 


dS/ 

Eq. 

(B-1) 


.1/5 


*The reference list is included at the end of this appendix; Table B-1 gives the 
nomenclature . 
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where according to reference (B-2) : 


(i®N 

V V, dS^ \y Pj )J 

oo ' — O — * 


1/2 


Eq. (B-2) 


The effective sweep angle A is given, in Reference B-3 , 


. -1 


Agjj - sin [sin A cos <» + sinacosAsini^] 


Eq. (B-3) 


with a being the angle of attack, the dihedral angle and A the actual sweep angle. 

For laminar flow with a freestream Mach number less than two, CAVE uses the 
modified Lees equation for the stagnation point: 


h = 0.5 
o 


Pr 


P 

2/3 




1/2 


Eq. (B-4) 


where the * quantities are evaluated at 


T* = 0.5 (T^ + T^) 


The velocity gradient term is computed using: 




1/2 


Eq. (B-5) 


For freestream Mach numbers greater than two, CAVE uses Detra, Kemp and 
Ridell's equation (Reference B-4) modified for effective sweep in accordance with 
References B-5 and B-6, to give the laminar heat transfer coefficient at the 
stagnation point: 


h_ = 


2.21 


o (T, 


\ 1000 / 


3. 15 


1.5 . 
cos A 


eff 


Eq. (B-6) 


116 



where A is the effective sweep angle given by Eq. (B-3). 

With the stagnation point coefficient h^ thus established, the distribution of 
convective coefficients around the leading edge is computed using Lees* formula 
(Reference B-7), with assumptions of cold wall and 


h = h F / 

O oo 


(dv^/dS)^ 


where 


F s 


V 


b f k ^'3 


- 1/2 


Maise and Rossi (Reference B-8) used this formula in subroutine LEES to obtain 
the distribution of h around the leading edge. The stagnation point value h^ was an 
unknown in their problem and was found by the CAPE code given the temperature 
history of the body. Basically, this same subroutine (LEESl) has been incorporated 
in the CAVE code. 

To evaluate F quantitatively, the local flow conditions around the leading edge 
are required. These conditions are found using the modified Newtonian law in the 
subsonic region surrounding the stagnation point and using Prandtl-Meyer expansions 
downstream of the sonic point locations. Specifically, in the subsonic region the 
local conditions will be predicted by the following equations : 


Po = Pt 


27 M^^sm‘0 - (7-1) 


2 

7+1 


( 7-1) sin^0 + 2 

' ' oo 


where 0 = 90“ - A 


eff 


The local pressure distribution is predicted according to Reference B-1: 


Pe " ^ (Pq ” Poo ) cos 0 


where 0 is the angle measured from the stagnation point. 


117 



“Pt 

*o ** 


(7+1) M^ sin^0 

7-1 

- 

7+1 

( 7 -1) M^ sin^+2 


2 7M^ sin^0 ~ ( 7 -1) 


T "r-l 


“e = 7 


^[ 0 ^ -31 


1/2 




/Ti^g 


e e 


P« = P^/RT 


e' e 


In the supersonic region downstream of the sonic points the local Mach number , 
Mg , is obtained via the Prandtl -Meyer turning angle v found from the following 
equations : 

■'2 = "l ^ [*2 - *l] 

where 


i/(M) 


= >/(7+l/Y-l) tan ^ V^(7 -1/7+1) (M^^-1) -tan”^ >/ m^-1 


and 



represents the change in flow angle . 
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With the local Mach number determined, the other pertinent flow properties are 
established using the following equations : 


P 


e 



(1 + 


y-1 

2 



-7/7-1 


o 


e V ^ e 


V = M a 
e e e 


= p / RT 
e e 


Flat Plate Geometry 


/ 
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For flat plate geometries Eckert^s reference enthalpy method (Reference B-8) 
is employed to predict the heat transfer rate , i.e. , 


h =. 


0.332p*V C 


e P 




(laminar flow) 


Eq. (B-7) 


and 


0.0296 p * V C 
h= — 


(turbulent flow) 


Eq. (B-8) 


The Reynolds number is based on the boundary layer length Jl . 

The subroutine FLATH calculates the local flow values to be used in the above 
equations based on oblique shock wave theory. That is, the shock angle 0 is found from 
solving the following cubic equation: 


7? + bZ^ + cZ + d = 0 


where , 


« 2 
Z = sin 9 

+2 



u - « 

A, • * 

- 7 Sin a 

c 

8 

2M^ + 1 



00 

C = : 

+ 

(7+X)' 


M 


r-1 


M 


CO_| 


. 2^ 
sin or 


d = - 


2 

cos a 


M 


oo 
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With the appropriate 9 selected for a weak shock, the other flow properties are 
calculated via the following equations : 




[- 


27 m 1 sin^a - (7-1) 


7+1 


] 


Tg = T, 


= V„ 

6 


[2 7 sin^a - (7-1)] [(7-1) sin^a + 2] 




(7+1)^ sin^a 


4 (M^ ^ sin^a -1) (7M^ sin^a + 1) 

(7+1)^ sin^a 


1/2 


Pe = Pe 


The subroutine TRANS establishes the laminar and turbulent regions using the 
boundary layer transition criterion shown in Figure B-1 (Reference B-9) . It is recog- 
nized that this criterion is not the final word on simple boundary layer transition 
criteria but rather representative of the best presently available . In light of this , 
the code is written so that as newer criteria are evolved they may be readily incor- 
porated into subroutine TRANS. For that matter, the entire aerodynamic heating 
portion of the code is written so that as the need arises to use more specialized 
equations , they may be easily substituted for those presented herein. 
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TABLE B-1. NOMENCLATURE AND UNITS 


Symbol 


Units 

a 

Speed of sound 

ft/sec 

<=P 

Specific heat at constant pressure 

Btu/lbm-^R 

d 

Hydraulic diameter 

ft 

9c 

Newton constant 

32.17 ft-lbm/lbf-sec^ 

h 

Convective heat transfer coefficient 

Btu/ft^-sec-°R 

k 

Thermal conductivity 

Btu-ft/ft^-sec*°R 

£ 

Effective boundary layer length 

ft 

M 

Mach number 

dimensionless 

P 

Pressure 

Ibf/ft^ 

Pr 

Prandtl number 

dimensionless 

q 

Heat flux 

Btu/ft^-sec 

r 

Recovery factor 

dimensionless 

R 

Gas constant 

lbf-ft/lbm-°R 

Re 

Reynolds number 

dimensionless 

Rfyj 

Nose radius 

ft 

S 

Distance along surface 

ft 

T 

Temperature 

°R 

V 

Velocity 

ft/sec 

a 

Angle of attack 

degrees 

7 

Ratio of specific heats 

dimensionless 

0 

Angular position from stagnation point or oblique shock angle 

degrees 

A 

Wing sweep angle 

degrees 

P 

Mass density 

Ibm/cu-ft 

P 

Viscosity 

Ibm/ft-hr 

Subscripts 



AW 

Adiabatic wail 


e 

Property at edge of boundary layer 


eff 

Effective 



(Continued) 
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TABLE B-1. NOMENCLATURE AND UNITS (Cont'd) 


Symbol Units 

Subscripts (Cont'd) 

o Stagnation point 

r Recovery 

t Total condition (i.e., condition that would exist if fluid brought to rest isentropically) 

W Wail 

°° Free stream 

Superscripts 

* indicates fluid property evaluated at temperature T* given by 

T* = Tq + 0.5<Tyv - Tg) + 0.22 (T^ - T^) 
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APPENDIX C 


LINEARIZATION OF RADIATION COUPLING 


Leading edges, cooled panels and slab geometries may involve radiation heat 
transfer at the surface. Typically the direction of heat transfer is away from the 
body, serving to cool it. An exact treatment of radiation within each time step pre- 
cludes an eigenvalue -eigenvector solution which depends on the problem being a 
linear one . Thus radiation heat transfer is given a linear representation with each 
time step in such a way that for the step size approaching zero the exact solution 
is produced. This appendix presents the linearization within each time step. It is 
the common one of modifying the convection coupling to account for the radiation , 
which t 3 ^ically amounts to a reduction in the convection coupling . 

Consider a node diagram for a surface node of a body. 


where 


hA 

acA 

kA 

L 

T 

^AW 

T 

^R 


W 



= convection coupling 
= radiation coupling 

= conduction coupling 

= adiabatic wall temperature of fluid (more properly referred to as 
recovery temperature of fluid) 

= background radiation temperature (usually taken as 0°Kj 
= temperature of surface node 
= temperature of interior node 



The heat transfer into the surface due to convection and radiation is 


% = A - T^) . «A fr^ - T^) 

This may be rewritten as ; 

where 

h - ~ 

^AW-V 


Eq. (C-1) 


Eq. (C-2) 


Eq. (C-3) 


Eq. (C-2) can be rewritten into the form used with CAVE; 

= heff A - T^) Eq. (C-4) 

where 

hea Eq. (C-5) 

Since T„ is usually taken to be 0°K, h„ will be negative and, therefore, the 

K rt 

effective convective coupling will be less than the actual convective coupling h. 

The term h„ is a linearized radiative coupling and it is in essence a correction factor 
K. 

to the convection coupling. Frequently, this correction factor amounts to less than 
5% of the convection coupling and is, therefore, of no significant consequence. On the 
other hand it is possible , particularly for high altitude trajectories , to have a 
relatively large correction factor, so large that the effective convective coupling is 
actually negative. A similar situation can occur late in a re-entry vehicle trajectory 
following peak heating. A negative h does not appear to pose any particular difficulty 
to the matrix routines within CAVE. However, a difficulty has been observed when 
^eff zero (i.e. , radiation coupling and convection coupling essentially 

equal). In these situations there are ill-conditioned matrices involved in the prob- 
lem solution and the cumulative effect of arithmetic roundoff becomes a serious 
matter in the eigenvector-eigenvalue iteration procedure within subroutine IJEN. 
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This difficulty has been observed on occasion with the manifestation of a failure 
within IJEN to obtain estimates for the eigenvalues and eigenvectors that were within 
tolerance. (A statement to this effect is printed out by CAVE.) However, the temper- 
atures calculated appeared correct. This may be attributed to the eigenvector values 
being reasonable although not within tolerance and furthermore not being of great 
significance since the body is not changing very much in temperature due to the very 
small coupling between it and its environment. Double precision arithmetic would 
help to alleviate this ’’ill-conditioned” difficulty. 

One final item to be mentioned is that in calculating hj^from Eq. (C-3), CAVE 
uses the temperature of the surface nodes at the beginning of time subinterval (end 
of previous time subinterval) . Therefore , in problems where radiation plays an 
important role and where the surface temperatures are varying rapidly with time , the 
user should use small time subintervals (perhaps one fifth of the trajectory table 
time intervals). 
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APPENDIX D 


PROGRAMMER ORIENTED DOCUMENTATION OF THE CODE 

This appendix presents the details of the CAVE organization and structure. A 
simplified logic flow diagram of CAVE is given in Figure D-1. 

CAVE is organized in a main program with 36 subroutines. The list of sub- 
routines is given in Table D-1 together with the function of each subroutine and the 
calling routine. Figure D-2 presents the organization of CAVE in terms of the more 
important subroutine calls. 

Many subroutines are either identical to, or modified versions of, subroutines 
used in the CAPE code (NASA CR-2435), which was developed by Maise and Rossi to 
solve the Inverse problem of finding the convective heat transfer coefficient given 
temperature history information. A modified version of a CAPE subroutine has a 
slight change in name to avoid any possibility of the wrong subroutine being used, as 
would be the case if both the CAVE and CAPE codes were resident on the same disc. 
To illustrate typical name changes, consider the CAPE subroutines PCP, HETRA 
and SLAB; in CAVE these subroutines in modified form are referred to as PCP4, 
HETRAl and SLAB2, with no significance attributed to the integer value. 

Flow charts or descriptions for each subroutine are given in Figures D-3 
through D-26. For convenience and completeness, they are given for subroutines 
that are identical to the CAPE routines reported in NASA CR-2435 as well as the 
new subroutines. 

Sets of input data for check cases were presented in Sections 3 through 5 where 
the standard output for these were described. For detailed output that maps the 
iterations and eigenvalues , the flags LTE and MON must be set equal to +6 in state- 
ment cards within SIZE2 formal values are -6). The detailed output refers to the 
iterations performed in the subroutine DESDAl and the subroutines that are called 
from it. The detailed output can be used to establish whether sufficient eigenvalues 
have been selected by the user since the detailed output shows the contribution of 
each term in the equation (Refer to Appendix A, Eq. (A-8) with the number of 
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nodes N replaced by the number of dominant eigenvalues NE): 


NE 

= ^ S Cjjexp(Xjt) 

‘ 3=1 

Each succeeding contribution should be smaller than the previous one, with the con- 
tribution of the last term being small in comparison with the accuracy desired for the 
temperatures. If this is not the case, NE should be increased and the problem rerun 
at the expense of increased computer time. 

CAVE prints out self-explanatory diagnostic messages for some of the errors 
that may be caused by faulty input data preparation. A diagnostic message "FAILURE 
IN UEN" signifies a failure of the Jennings algorithm to converge within the maximum 
number of iterations in the subroutine UEN. This message has been observed on rare 
occasions as a result of a matrix being ill-conditioned at a particular time in the flight 
trajectory where the radiation -convection coupling is approaching zero (refer to 
Appendix C). In these situations the temperatures have appeared to be correct and 
the message has been ignored. If corrective action should be necessary it could in- 
clude the following: increasing the maximum number of iterations as given by NIJ 
in subroutine SIZE2 from its present value of 20; increasing the tolerance on the 
convergence test as given by TIJ in subroutine SIZE2 from its present value of 0.1; 
revising the time steps used to avoid this troublesome point in the flight trajectory; 
and as a last resort utilizing double precision arithmetic. 
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FIG. D-1 SIMPLIFIED FLOW OF LOGIC IN CAVE 
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Table D'1. Subroutines Used in CAVE 


Name 

Called By 

Main Purpose 

SIZE2* 

CAVE (main prog) 

Computes storage locations needed. Compares to number requested 

OVLYIO* 

SIZE2 

Sets up arrays for PCP4 

PCP4* 

OVLYIO 

Reads and writes property data, controls geometry 

SLAB2* 

PCP4 

Computes volumes and conduction shape factors for slab and cooled 
panel geometries 

LEAD4* 

PCP4 

Computes volumes and conduction shape factors for leading edge 
geometries 

X24C 

PCP4 

Computes volumes and conduction shape factors for basic X24C 
geometries 

GEN 

PCP4 

Reads and writes volumes and conduction shape factors for general 
geometries 

MATOUT 

PCP4 

Writes material properties 

XTABS1 

GEN 

Reads tabular values of hA and T^yy for general geometry problems 

OVLY20 

SIZE2 

Reads initial temperature distribution and flight trajectory. Controls 
problem solution, steps time, computes average convection couplings 
for each time step. Writes solution each time step 

LINFIT 

OVLY20,PROP 

Finds value from a table by linear interpolation 

ATTAC2* 

OVLY20 

Finds nude number closest to stagnation point and renumbers nodes 
as required by LEES1 

LEES1* 

ATTAC2 

Computes ratios h/h^p and variation around leading edge problems 

NURED 

OVLY20 

Reads tabular values of hA and as functions of distance and time 

DINTK 

OVLY20 

Finds value from a table using double interpolation 

PROP 

OVLY20 

Computes conduction couplings and mass specific heat product for 
each element given conduction shape factors and volumes 

PRPOUT 

OVLY20 

Writes node numbers, conductances, capacitances and initial 
temperatures 

XINTP1 

OVLY20 

Finds values of several dependent variables from a table by linear 
interpolation on a single independent variable 

FLATH 

OVLY20 

Finds h and T^y^ for flow over a flat plate 

ATMOS 

OVLY20,FLATH 

Finds atmospheric pressure, temperature and density for given altitude 

POLRT 

FLATH 

Computes the roots of a polynomial 


(continued) 
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Table D-1. Subroutines Used in CAVE (Cont'd) 


Name 

Called By 

Main Purpose 

TRANS 

FLATH 

Determines whether flow over flat plate will be considered laminar 
or turbulent 

DESDA1 

OVLY20 

Calls eigenvalue and matrix routines. Calculates temperatures 

IJEN+ 

DESDAl 

Obtains dominant eigenvectors and eigenvalues of a given matrix 
(using Jennings method of simultaneous vector iteration) 

EIGVC+ 

DESDA1 

Prepares approximate guesses for the eigenvectors to start the Jennings 
algorithm iteration for the first time step 

BFACS+ 

UEN 

Factorizes a banded positive-definite matrix 

BSOLS+ 

DESDAl, UEN 

Using the factors of a given banded positive-definite matrix A as 
generated by BFACS solves for X in the system AX = Y 

ORNML+ 

UEN, DESDAl 

Carries out the standard Gram-Schmidt orthonormalization of a 
group of vectors 

HETRA1* 

DESDAl 

Sets up coefficient matrix (of conductances) in compact form 

RVORDR+ 

UEN 

Reorders estimated eigenvalues according to magnitude 

AORDER+ 

UEN, RVORDER 

Sets up permutation indices needed for ordering the eigenvalues 

DISPLA+ 

(Various) 

Prints information, mainly debug special output, in array form 

PART+ 

DESDAl, PCP4 

Prints debugoutput information and intermediate timing of calculation 

DATE+ 

PART 

Determines data of run 

SWITCH+ 

DISPLAY 

Converts columns of a matrix to rows or visa versa 

SCAPR2+ 

(Various) 

Computes scalar product of two vectors 


*Modified version of CAPE routine + CAPE routine 
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CAVE: ASSIGN TOTAL 
DIMENSION 




GEN: (GENERAL GEOMETRY) 

READ: VOLUMES 

CONDUCTION SHAPE FACTORS 
CONVECTION COUPLINGS 
ADIABATIC WALL T 


0 


SLAB2: (SLAB & COOLED PANEL 
GEOMETRIES) 

READ & WRITE AX's, AY's TAU, SI , ETC. 
COMPUTE CONDUCTION SHAPE FACTORS 
COMPUTE VOLUMES 
COMPUTE CONVECTION COUPLINGS 
FOR COOLED SURFACES 


X24C: (X24C STRUCT! 

READ & WRITE AX's, A 
COMPUTE CONDUCTIO 
COMPUTE VOLUMES 
COMPUTE CONVECTIOr 
COOLED SURFACES 

JRAL ARRANGEMENT) 
V's,S1,S2, ETC. 

M SHAPE FACTORS 

M COUPLINGS FOR 

1 


^ 

r 

LEAD4: (LEADING EDGE GEOMETRY) 

READ & WRITE MCAP, THETA, AX's, AY's, ETC, 
COMPUTE CONDUCTION SHAPE FACTORS 
COMPUTE VOLUMES 

COMPUTE CONVECTION COUPLINGS FOR 
COOLED SURFACES 


FIG. D-2 ORGANIZATION OF CAVE IN TERMS OF THE MORE IMPORTANT SUBROUTINE CALLS 
(SHEET 1 OF 2) 
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OVLY20; 




READ INITIAL T DISTRIBUTION 
READ TRAJECTORY INPUTS 
READ TIME INTERVALS 
READ HEATING MODIFIERS 





FIRST TIME STEP 
OR 

T DEPENDENT PROPERTY? 



< 


FIRST TIME STEP? 



SLAB, COOLED PANEL 
OR X24C GEOMETRY 
WITH h COMPUTED 
BY CAVE? 



LEADING EDGE GEOMETRY 
WITH h COMPUTED BY 
CAVE? 




FIND TEMPERATURES 
AT END OF THIS 
TIME INTERVAL 


WRITE: 


h 


hA 


Taw 


T 


Tss 



TIME EQUAL 
FINAL TIME? 


YES 


YES 


PROP: 


COMPUTE THERMAL CAPACITANCE 
AND CONDUCTION COUPLINGS 


=r 


YES 


PRPOUT; WRITE NODE NUMBERS, THERMAL 
CAPACITANCES, CONDUCTION COUPLINGS 


YES 


FLATH: COMPUTE CONVECTIVE 

COEFFICIENT AND ADI ABATIC WALL 
TEMPERATURE FOR FLAT PLATE SURFACE 


YES 


ATTAC2: FIND NODE 
NUMBER CLOSEST TO 
STAGNATION POINT AND 
RENUMBER NODES FOR 
LEES1 


LEES1: COMPUTE 
h/hsp& Taw 
AROUND LEADING 
EDGE 


6 


DESDA1: 

• FIND EIGENVALUES 
AND EIGENVECTORS 


• COMPUTE STEADY- 
STATE TEMPERA- 
TURES 

• COMPUTE TEMPER- 
ATURES END OF 
THIS TIME 
INTERVAL 


HETRA1: SETS UP 
COEFFICIENT MATRIX 
IN COMPACT FORM 


IJEN: FINDS EIGEN- 

VALUES & EIGENVECTORS! 


RVORDR: REORDER 
EIGENVALUES ACCORD- 
ING TO MAGNITUDE 


BSOLS; SOLVES 
MATRIX EQUATION 
AX = Y 


SCAPR2: COMPUTES 

SCALAR VECTOR 
PRODUCT 


FIG. D-2 ORGANIZATION OF CAVE IN TERMS OF THE MORE IMPORTANT SUBROUTINE CALLS 
(SHEET 2 OF 2 ) 
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FIG. D-3 SUBROUTINE SIZE2 FLOW CHART 
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FIG. SUBROUTINE PCP4 FLOW CHART 










FiG. D-5 SUBROUTINE SLAB2 FLOW CHART 


139 










FIG. SUBROUTINE LEAD4 FLOW CHART 
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FIG. D~7 SUBROUTINE X24C FLOW CHART 





GEN 
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GEN FLOW CHART 









FIG. D-9 SUBROUTINE MATOUT FLOW CHART 





















FIG. D-11 SUBROUTINE OVLY20 FLOW CHART (SHEET 2 OF 4) 
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FIG. D-1 1 SUBROUTINE OVLY20 FLOW CHART ISHEET 3 OF 4) 
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FIG. D-11 SUBROUTINE OVLY20 FLOW CHART (SHEET 4 OF 4) 
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FIG. D-13 SUBROUTINE ATTAC2 FLOW CHART 
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f IG. D-15 SUBROUTINE PROP FLOW CHART 




ENTER 


THERMAL^ 

CAPACITANCE 

EQUAL 

V ZERO ^ 









FIG. D-17 SUBROUTINE XINTP1 FLOW CHART 


154 




























I 7'i' 



SET TURBULENT 
FLOW INDICATOR 


SET LAMINAR 

. 

FLOW INDICATOR 




►{RETURN 




COMPUTE 
TRANSITION 
REYNOLDS 
NUMBER, RET 



COMPUTE 

NO 

TRANSITION 


REYNOLDS 


NUMBER, RET 


COMPUTE 
TRANSITION 
REYNOLDS 
NUMBER, RET 


YES 

SET 

TURBULENT 


FLOW 


INDICATOR 


1 

r 




FIG. D-20 SUBROUTINE TRANS FLOW CHART 


15T 
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FIG. D-21 SUBROUTINE DESDA1 FLOWCHART 
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FIG. 0*22 SUBROUTINE IJEN FLOW CHART 












•(SOURCE: NASA CR‘2435) 


FIG. D-23 SUBROUTINE ORNML FLOW CHART* 
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FIG. D-24 SUBROUTINE HETRAI FL6W CHART 




QUANTITY 

SYMBOL 

INPUT/OUTPUT 

DIMENSION 

EIGENVALUES 

R 

IN + OUT 

R(MM) 

EIGENVECTORS 

V 

IN + OUT 

V(MID,MM) 

RANK VECTOR 

K 

OUT 

K(MM) 

PERMUTATION VECTOR 

L 

OUT 

L(MM) 

EIGENVECTOR DIMENSION 

N 

IN 

— 

NUMBER OF EIGENVECTORS 

MM 

IN 


DIMENSION OF ARRAY 

USED TO STORE EIGENVECTORS 

MID 

IN 




•(SOURCE: NASA CR-2435) 

FIG. D-25 SUBROUTINE RVORDR FLOW CHART* 
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III ll■l■l■ll ll■■lll 






SCAPR2 CALCULATES THE INNER 

PRODUCT OF TWO VECTORS STORED 

AS EQUALLY SPACED WORDS IN FORTRAN ARRAYS 



•(SOURCE: NASA CR-2435) 


FIG. 0-26 SUBROUTINE SCAPR2 FLOW CHART* 
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Purpose 


Sheet D'1. Subroutine NURED Description (Sheet 1 of 4) 


To read a set of tables for functions 2 variables 


Program Description 

The input data must be structured as specified below. The calling program must contain the statement, 
COMMON STG(L), where L >{L1(1)+1) {L2(1)+1) + ... + (L1(M)+1) (L2(M)+1). Program restrictions 
are noted below: 


Input Parameters 

FORTRAN Name 


Description 


NUMTBL 


MANDAN 

LI 


L2 

NUMPTS 


Output Parameters 
NG 


Calling Sequence 


= 1 for first call to NURED 

= K for replacement of tables, where K is the table number of the first table 
being replaced 

= 0 for initial read in 
= 1 for table replacement 

Array of dimension M, where M is the maximum number of tables in storage 
at any one time 

LI ( K) = number of x j in table K 

Array of dimension M 

L2(K) = number of yj in table K 

Array of dimension M+1 

NUMPTS(K) = the number of table entries preceding table K 


To be set to 0 in the calling program 
NG = 0 if tables have been read correctly 
NG = 1 if there has been a read in error 


CALL NURED (NUMTBL, MANDAN, NG, LI, L2, NUMPTS) 


Input Format 

For each table the functions values, F, (X,Y) are entered for consecutive groups of nine values of X, over all 
values of Y 


(continued) 
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Sheet D-1. Subroutine NURED Description (Sheet 2 of 4) 


Input Format (Cont'd) 

The input data must be structured as indicated. First, there is a header card with the appropriate entries in 


the indicated columns 





1-2 

3-4 



17 

70 

71-72 

LI 

L2 



COMMENTS 


Seq # 

Next, for the first 9 values of X the input data takes the form 



1-7 

8-14 

15-21 

22-28 

. . . 57-63 

64-70 

71-72 


Xl 

X2 


... Xg 

Xg 

Seq # 

Yl 

F(1.1) 

F(2,1) 

F(3.1) 

... F(8,1) 

F(9,1) 

Seq # 

^^L2 

F(1,L2) 

F(2,L2) 

F(3,L2) 

. . . F(8,L2) 

F(9,L2) 

Seq# 

For the second group 

of 9 values of X, the form of the input is 




^10 

^11 

Xi2 

... X^7 

Xl8 

Seq # 


F(10,1) 

F(11,1) 

F(12,1) 

... F(17,1) 

F(18,1) 

Seq # 

^L2 

F(10,L2) 

F(11,L2) 

F(12,L2) 

... F(17,L2) 

F(18,L2) 

Seq# 


Additional values of X are handled similarly. When all X values have been accounted for, the next table 
of values for the next function is set up in the same way. When ail data tables which are to be read in at 
any one time have been set up as indicated, a blank card is placed at the end of the data deck. 

The parameters used above are defined as follows: 


LI = Number of X values in table (12) 

L2 = Number of Y values in table (12) 

Seq # = Sequence number of a card within a table, beginning with 0 for the first' card (12) 

F(I,J) = Function value for X|, Yj (E7.0) 

X| = Argument 1 values in table (E7.0) 

Y| = Argument 2 values In table (E7.0) 


The figure illustrates the card format for the tables (see Sheet 4 of Sheet D-1 ). 


(continued) 
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Sheet D-1. Subroutine NURED Description (Sheet 3 of 4} 


Program Restrictions 

The tabular values of the variables X and Y must appear in algebraically increasing order. The variables 
X and Y, and the function values must: 

• Be single precision numbers less than 99999E9 

• Have a maximum of 7 significant digits if positive 

• Have a maximum of 6 significant digits if negative 

A maximum of 99 cards is allowed for each table 
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CONTINUATION CARDS FOR REMAINING Y VALUES 



SHEET D-1 SUBROUTINE NURED DESCRIPTION (SHEET 4 of 4) 
















































Sheet D-2. Subroutine DINTK Description (Sheet 1 of 3) 

Purpose 

Table lookup and linear interpolation for several functions of two variables. 

Analytic Description 

In the derivation of STINT (6.1. 1.5) equations, we find the linear interpolation form 

f(x,1) Iv - yn) + f|x,0) (vi - y) . uu i u • i ♦ 

yl _ y t which, by algebraic manipulation becomes 

Vi ’ Vo 


(i) 


f(x.y) 


(y^ - y) (yo“ y) 

vi-vo 


Analytic Restrictions 


f(x,0) f(x,1) 

Vq-V Vt-VL 


The function f should be linear 


Program Description 

f-| (x,y), f2<x,y), ..., f^ (x,y) are found for Xj < x < Xj+^ and Vj < y < yj+^ , using equation (i) 

Program Restrictions 

Extrapolation will not be performed. The calling program must contain the statement COMMON STG(L>, 
where 


L > 


2 

:=1 


(Lid) + 1) (L2(l) +1) 


Input Parameters 

FORTRAN Name 


Description 


LI, L2, NUMPTS 

KODE 

N1HIB4 

N2HIB4 

ARG1 

ARG2 

NUMTBL 

L3 


As described in NURED 

Dummy array of length m. Initialize to zero in the calling program 

Dummy array of length m 

Dummy array of length m 

Value of x argument 

Value of y argument 

Number of the first table in which interpolation will take place 
Number of functions to be interpolated 

(continued) 
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Sheet D-2. Subroutine DiNTK Description (Sheet 2 of 3) 

Output Parameters 

FCT Arrayof length L3, consisting of the elements FCT(1) = Fj (x,y), 

FCT(L3) = f„(x,y). (pNUMTBL, n=i+L3-l) 

NG To be initialized to zero before calling DINTK 

NG = 0 indicates a normal return 

NG = 2 indicates a machine error or an error in the calling sequence 
NG = 3 indicates that x or y is outside of the range of the tables 

Library Supplied Routines 

User must call NURED to read in the tables of x, y, f-|(x,y), f^(x,y) prior to calling DINTK 
Calling Sequence 

CALL DINTK (LI, L2. NUMPTS, KODE, N1HIB4, N2HIB4, ARG1, ARG2, NUMTBL, L3, FCT, NG) 


List of Variables 


FORTRAN Name Description 


LI 

L2 

NUMPTS 

KODE 

N1HIB4 

N2HIB4 

ARG1 

ARG2 

NUMTBL 

L3 

FCT 

NG 

STG 

X 

Y 

ANS 

FACTOR 

•= 

NT 

NCHECK1 

NCHECK2 

INC 

K 

NL01 

NL02 

NBL01 

NBL02 


Table of values of input 
Temporary storage for ARG1 
Temporary storage for ARG2 
Working storage 
Working storage 
Working storage 


Indices 


(continued) 


Sheet D-2. Subroutine DINTK Description (Cont'd) (Sheet 3 of 3) 


List of Variables (Cont'd) 


FORTRAN Name 


NHI1 

NH(2 

JINDEX 

KINDEX 

I 

J 

M 

LDUMY1 

LDUMY2 

LENGTH 

LIMLOW 

LIMUPR 

KTEST1 

KTEST2 

KTEST3 

KTEST4 

L 

M 

I 

J 

LA 

LB 

NDUMY1 

NDUMY2 




> 


J 




> 


J 



Description 


indices 


DO loop indices 


Working storage 
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Sheet D-3. Subroutine POLRT Description 


Purpose 

Computes the real and complex roots of a real polynominal 


Usage 


CallPOLRT(XCOF.COF.M.ROOTR,ROOTI,IER) 

Description of Parameters 

XCOF = Vector of M+1 coefficients of the polynominal ordered from smallest to largest power 
COF = Working vector of length M+1 
M = Order of polynominal 

ROOTR = Resultant vector of length M containing real roots of the polynominal 

ROOTI = Resultant vector of length M containing the corresponding imaginary roots of the polynominal 

lER = Error code where: 

IER=0 No error 
I ER=1 M less than one 

IER=2 M greater than 36 

IER=3 Unable to determine root with 500 interations on 5 starting values 
IER=4 High order coefficient is zero 

Remarks 

Limited to 36th order polynominal or less 

Floating point overflow may occur for high order polynominals but will not affect the accuracy of the results 
Subroutines and Function^ubprograms Required 
None 
Method 

Newton-Raphson iterative technique. The final iterations on each root are performed using the original 
polynominal rather than the reduced polynominal to avoid accumulated errors in the reduced polynominal 
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Sheet D-4. Subroutine EIGVC Description (From NASA CR-2435) 

EIGVC computes guesses of the eigenvalues, eigenvectors and associated permutation index that are necessary to 
start the iteration in the Jennings method to calculate eigenvalues and eigenvectors. The formulae used for these 
guesses are: 


ith eigenvalue R - i 

ith eigenvector A ^ e ■’'V/{h/2) sin| + 7r-~ (i - T) | 


y'T” FOR TOP 
,, HALF 


f FOR BOTTOM 
y I HALF 
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Sheet D-5. Subroutine BFACS Description (Sheet 1 of 2) (From NASA CR-2435) 


Given a special, L-banded, positive or negative definite symmetric N-th order matrix, S, decompose it into the 
product: 

S = U ' D~* U 


where U is an L-banded nonsingular upper triangular matrix with unit diagonal elements, and D is a nonsingular 
diagonal matrix. 

The S matrix is inputted as elements of three one-dimensional arrays. A, B, E. The N elements of A are the main 
diagonal elements of S, the leading (N— L) elements of E are the L-th super- (and sub-) diagonal elements of S. The 
N-1 elements of B are super- (and sub-) diagonal elements of S. *The trailing L/2 elements of E (optionally), define 
main cross diagonal elements of the upper L-th order submatrix of S. In general, later definitions override earlier one, 
e.g. if L = 1, B (not E) defines the super diagonal elements of S. In the case for which BFACS is intended, most elements 
inside the band are zero. The cross diagonal is installed only if the argument, a, is not zero. The S matrix is topolo- 
gically equivalent to conduction paths in a slab (leading edge if a 0); consider the N = 12, L = 4 example: 


1 ■ 

5 

• 9 

2 

6 

■10 

3 ■ 

■ 7 

•l'l 

4 • 

■ 8 ■ 

■12 


Note: b^ = bg = 0 
for the conduction 
problem generally 
but must be explicitly 
made 0 for BFACS 

Note: If nose paths 
are included 

are used (q 0). 



1 

2 

3 

4 

\ 

5 

1 

6 

7 

8 

9 

10 

11 

12 

1 

®1 

b^ 0 

1 ®12 
1 

1 ®1 
1 

0 

0 

0 

0 

0 

0 

0 

2 

bl 

^2 1 
1 1 

®1 1 ' 
1 

' 0 
1 


0 

0 

0 

0 

0 

0 

3 

0 

1 ®11 ' 
v" 

j 

0 

^3 

^3 

1° 

0 

®3 

0 

0 

0 

0 

0 

4 

®12 

^3 

^4 


0 

0 

®4 

0 

0 

0 

0 

5 

«1 

0 

® 1 "5 

^5 

0 

_ 

“ 1 

®5 

0 

0 

0 

6 

0 

®2 

0 

0 

>^5 

1 

^6 

^6 

0 1 
1 

0 

®6 

0 

0 

7 

0 

0 

®3 

0 

1 

i“ 

^6 

^7 

"7 1 

0 

0 

®7 

0 

8 

0 

0 

0 

®4 

1 0 

1 

0 

h 

— r 
00 
(D 

1 


0 

0 

®8 

9 

0 

0 

0 

0 

® 

0 

0 

®\ 

^9 

^9 

0 

0 

10 

0 

0 

□ 

0 

0 

® 

0 

0 1 

^9 

^10 

^10 

0 

11 

0 

0 

0 

0 

0 

0 

«7 

0 

0 

^10 

®11 

bll 

12 

0 

0 

0 

0 

0 

0 

0 

eg 1 

0 

0, 

bl1 

ai2 


173 




Sheet D-5. Subroutine BFACS Description (Sheet 2 of 2) (From NASA CR-2435) 


To take advantage of both symmetry and the band form U is stored in a rectangular array of size MID by L, 
where MID > N; the bottom row is used as scratch storage, r,|, r2< . • • , r|_, and the unused bottom triangle is 
zeroed out (for convenience in printing only). U appears as: 

gyQP^Q array non zero elements of u matrix 


234 123456789 10 n 12 



_ 


- " 

1 

^11 ^13 ^14 

1 

1 Ui2 u^4 

2 

“21 “22 “23 “24 

2 

1 U21 U22 ^23 U24 

3 

“31 “32 “33 “34 

3 

1 U3^ U32 U33 U34 

4 

“41 “42 “43 “44 

4 

1 U4^ U42 U43 U44 

5 

“51 “52 “53 “54 

5 

^ “51 “52 “53 “54 

6 

“61 ^62 “63 “64 

6 

1 Ugi Ug2 Ug3 Ug4 

7 

“71 “72 “73 “74 

7 

1 U^^ U72 U73 U74 

8 

“81 “82 “83 “84 

8 

1 Ugi Ug2 ug3 Ug4 

9 

“91 “92 “93 ® 

9 

1 Ug^ Uq 2 Ug3 

10 

“10,1 “10,2® ® 

10 

^ “10,1 “102 

11 

u^ 1 0 0 0 

11 

^ “11,1 

12 

ri rj rg 

__ — 

12 

1 


The N elements of D"^ are stored in an N-array. For the usual case of S being either positive-or negative-definite, 
these elements are all positive or all negative, respectively. However the routine will "work” provided only that the 
leading N principle minors are non-zero. For details see the following article which guarantees high accuracy only for 
the definite cases of usual interest; "Symmetric Decomposition of Positive Definite Band Matrices R.S. Martin, 

J.H. Wilkinson, C. 1/4, LINEAR ALGEBRA - HANDBOOK FOR AUTOMATIC COMPUTATION, VOLUME II, 
Springer-Verlag, 1971. 
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Sheet D-6. Subroutine BSOLS Description (From NASA CR-2435) 


Given the product form decomposition of an L-banded symmetric matrix, S = U^D*’ U, as calculated by the 
BFACS subroutine, BSOLS solves a system of N linear equations with M right hand sides: 

S jY, Y2 - YM| = j Y, Y2 - • YMj 

The routine si mply carries out the standard forward substitution phase: 
z = U"T y 

followed by the standard backward substitution phase: 

X = U' D z 

The only unusual aspect is the rather unorthodox storage scheme which is described in the documentation for 
subroutine BFACS. This scheme is necessary to exploit the banded symmetric form of S in the most efficient way 
in terms of computer memory. For details see: R.S. Martin, J.H. Wilkinson, "Symmetric Decomposition of Positive 
Definite Bank Matrices", in: Linear Algebra— Handbook for Automatic Computation, Volume M, C. 1/4, Springer- 
Verlug, 1971 
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Sheet D-7. Subroutine AORDER Description (From NASA CR-2435) 


PURPOSE: ORDER A SET OF REAL NUMBERS 


CALLING SEQUENCE: CALL AORDER (A, N, IPERM) 


NAME DIMENSION DESCRIPTION 


INPUT A 

N 


A(N) ELEMENTS TO BE ORDERED 

NUMBER OF ELEMENTS = INI 
N >0 INCREASING ORDER 
N < 0 DECREASING ORDER 


OUTPUT IPERM IPERM(N) 


ORDER VECTOR - 
SPECIFIES THE SEQUENCE 
OF ELEMENT INDEX NUMBERS 
WHICH WILL PRESENT 
A AS AN ORDERED SET, 
i.e. 

DO 100 I = 1, N 
100 WRITE (6,1) A (IPERM(D) 

1 FORMAT (F 10.5) 


WILL LIST A AS AN ORDERED ARRAY 


AORDER CALLS NO OTHER SUBROUTINES 



Sheet D-8. Subroutine DISPLA Description (From NASA CR-2435) 


TITLE: DISPLA - Prints scalars, vectors, rectangular matrices, packed symmetric matrices, and Hessenberg 

matrices. 

AUTHOR: M.J. Rossi 

DATE: September 1973 

APPLICABLE COMPUTERS; IBM 360/370; CDC 6000 SERI ES 

SOURCE LANGUAGE: FORTRAN IV 


PURPOSE: To simplify printing of mathematical types of data structures in an easily read formal which allows 

titles and index labels. 


METHOD: FORTRAN looping and write statements which indexes and addresses arrays according to their type. 

USAGE: Call DISPLA {X, NFILE, TITLE. KAR, KIND, NROWS. NCOLS, MID). 


X — Input — Array of one or more values to be printed 
NFI LE - Input — FORTRAN unit for printing. 

TITLE — Input — Vector of KAR characters used as title. 

KAR — Input — Number of characters in above string. 

KIND — Input — Type of mathematical data structure: 

= O scalar (or vector printed on one line with no index) 

= 1 vector of iNROVlfSI elements, indexed 

= 2 Rectangular INROWSI by NCOLS matrix — Dimension (MID, *) 
= 3 Packed Symmetric matrix of order | NROWS | 


1 

2 

4 


2 4 

3 5 
5 6 


— lower triangular partial rows if NROWS positive 


1 

2 

3 


2 

4 

5 


3 

5 

6 


— lower triangular partial columns if NROWS negative 


= 4 — Transposed Hessenberg matrix of order NROWS — Dimension (MID, MID) 


NROWS 

— Input 

NCOLS 

— Input 

MID 

- Input 


Number of elements if KIND = 0 or 1 

— Number of rows if KIND = 2 

— Matrix order if KIND = 3 or 4 

— Number of columns if KIND = 2 

— Ignored-otherwlse 

— Matrix Dimension if KIND = 2 or 4 

— Ignored otherwise 


SUBROUTINE REQUIRED: SWITCH 
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Shest D-9. Subroutine PART Description (Prom NASA CR-2435) 

TITLE: PART — Prints standard 120 character labels at the top of the next page. 

AUTHOR: M.J. Rossi 

DATE: September 1973 

APPLICABLE COMPUTERS: IBM 360/370; CDC 6000 SERIES 

SOURCE LANGUAGE: FORTRAN IV with 2 Assembler Language Subordinate Subroutines. 

To make it convenient to produce standard printed labels with "part" numbers, date, and running CPU 
time on the top of the next page. Also, prints short line on next line with just CPU time for intermediate 
timing. • , 

On the first printing entry for a given computer run the Date Subroutine is invoked and an 8-character 
field of an internal word is stored with the date in the form: "KK/LL/MM", where KK is the index 
number for the month, ,LL is the day of the month, and MM is the last 2 digits of the year, e.g., 

March 15, 1973 3/15/73. Also, at this time, the SECOND subroutine is invoked to both 

establisft the zero time point and to set the units to hundreds of a second. Then the first printed page 
heading is given with a zero time and a PART number 1 reported. Subsequent printing calls will give 
the time as; NN.ILJJ where NN is the number of minutes elapsed, II is the number of seconds, and 
JJ is the number of hundredths of seconds. The PART number is incremented by one for each 
printing call. There are two fields of alphameric information for the full printing mode which are under 
control of the user: (1} The first is a 40 character LABEL field which is set upon calling PART in the 
non-printing mode. \ 2 ) The second is a 48 character field which is supplied on a full printing cal). 

There is also a partial printing mode which simply results in the appearance on the next line of an 8 
character field of user supplied TITLE along with running CPU time. 

Calf PART I'XX. . -X'. I 1 

. ..X' — Input = Alphameric string of either 8, 40, or 48 characters depending on the value of L. 
t FORTRAN unit for printing, if positive 

' • — If zero, simply sets 40 character LABEL field and returns 

— If negative, prints 8 character TITLE — 'XX X' — and CPU time on next line and 

increments PART number. 

— If positive, prints DATE, TIME, 40 character LABEL, 48 character TITLE, Part Number 
and spacers with standard notation. 

SUBROfftC^ESJ^TlEQUIRED; DATE, SECOND 


PURPOSE: 

METHOD: 

USAGE: 
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Sheet D-10 Subroutine Switch Description (From NASA CR-2435) 


PROGRAM TITLE: Utility routine for re-arrangement of certain triangular arrays 


SUBROUTINE NAME: SWITCH INDEX: 

ANALYST; F. Nolan 

PROGRAMMER: F. Nolan DATE; 

DOCUMENTATION AUTHOR; F. Nolan DATE; 

SOURCE LANGUAGE; FORTRAN IV 

APPLICABLE COMPUTERSr IBM Systems 360, 370;CDC 6000series 
REVISION; DATE: 


12.6.0.1 

June 15, 1967 
June 20. 1973 


PURPOSE; To provide a convenient conversion between two common arrangements for the storage of triangular 
(and symmetric) matrices. 


ANALYTIC DESCRIPTION: The routine makes systematic use of transpositions, i.e., interchanges of two array 

elements. It is a well known result in permutation theory that every permutation can be represented 
as a product of transpositions. 


PROGRAM DESCRIPTION: There is no loss of generality in assuming that the input matrix is of lower triangular 

form. It is natural to store such matrices by row or by column. Both arrangements are illustrated for 
a matrix of order 5. The understanding is that the (4,3) element, for example, is assigned position 9 
using row storage, and position 1 1 using column storage. 


Row Storage 

Column Storage 


1 

1 




2 3 

2 

6 



4 5 6 

3 

7 

10 


7 8 9 10 

4 

8 

11 

13 

11 12 13 14 15 

5 

9 

12 

14 15 


Given a lower triangular or symmetric matrix, stored in either fashion, SWITCH can re-arrange it to the 
other form. The re-arrangement is carried out "in place" in the sense that no auxiliary array is required. 
For an input matrix of order m, the transformation is performed in approximately transpositions. 
There are no rounding errors. 

PROGRAM RESTRICTIONS: The matrix must be of order at least 3. 

INPUT PARAMETERS: 

FORTRAN Name Description 


A Singly-dimensioned real array containing the matrix to be re-arranged. 


M 


Order of matrix A. If M is given positive, conversion is from row to column 
storage. If M is given negative, conversion is from column to row storage. 


OUTPUT PARAMETERS: 


A 


Matrix in re-arranged order. 


CALLING SEQUENCE: 


CALL SWITCH (A, M) 
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APPENDIX E 


DISCUSSION OF NONLINE ARITIES AND TIME DEPENDENCY OF h AND T * 

AW 

The heat conduction problem that CAVE solves generally contains nonlineari- 
ties as a consequence of the material properties being temperature dependent and the 
radiation heat transfer taking place at the surface of the structure . The problem 
usually has the further complication of time varying boundary conditions by virtue of 
h and being time dependent . This appendix considers two questions ; (1) how 
should the nonlinearities be handled; and (2) how should the time dependence of h and 

Taw be handled. 

AW 

Considering the second question first, it is noted that the time dependence does 
not destroy the linearity of the problem. However, in this case an eigenvalue solution 
cannot be obtained by treating h(t) directly in the system of ordinary differential 
-equations . This is seen by considering the restrictive case in which 

h(y.t) =\(y) (t) 

If the semidiscrete problem has the following form (Refer to eq. A- la) 

= A(t) T+b(t): T(0) = T^ 

where for simplicity M is taken to be the Identity matrix , then the diagonal elements 
of A depend on t and the above assumption on the form of h(y ,t ) leads to 

A(t) = A^-Ha^(t)J 

where £L^(t) is a scalar time function and J is a diagonal matrix with ones and zeros. 
Ih this case 

A(y A(t2> fi A(t2> A(tj) 

which therefore precludes an eigenvalue/eigenfunction solution. 


♦This appendix courtesy of Michael J. Rossi. 
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Iherefore, since the eigenvalue solution does not exist for even this restrictive 
form on h(y,t), the recommendation then Is to subdivide the time interval and take h 
and to be piecewise constant within each subinterval. Returning to the first of 
the two questions posed, we see another motive for subdividing the time interval. If 
the subdivision for h and produces short enough intervals one may account for 
the temperature -dependent properties and radiation terms by taking them to be piece- 
wise constant depending on the temperature at the beginning of each subinterval. One 
must caution, however, against taking time subintervals which are so short as to 
require too large a number of eigenvalues and eigenfunctions in order to ensure an 
accurate solution of the resulting subproblems. The best approach might be to first 
predict the temperature history on a subinterval based on a time-invariant linear 
model, and then. to correct the solution by considering a forcing term to accoimt for 
the total neglected effects of radiation, temperature dependent thermal properties, 
and time-dependent convection. 
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APPENDIX F 


DERIVATION OF SOLUTION TO THE EQUATION MT - BT +F 
This appendix presents a step-by-step solution to the vector equation 
MT = BT + F Eq. (F-1) 


which represents the temperature response within a body that has been discretized 
into a number of uniform temperature systems or nodes that are coupled and changing 
in temperature. The elements of the diagonal matrix M represent the mass-specific 
heat product of each system or node; the elements of the matrix B represent the 
convective -conductive couplings between nodes; and the elements of the vector F 
represent the product of the convective couplings with the corresponding fluid 
adiabatic wall temperature (or recovery temperature). And, of course, the elements 
of the vector T represent the instantaneous temperatures of the nodes ; the 
elements of the vector T the time rate of change of the temperatures. Specifically, 
for the four node system shown in Figure F-1: 


M « 


0 

0 


0>vCp)^ 

0 

0 


(pVCp>3 


(pVCp)^ 


1 = < 
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B = 


-(H+K^g) 

Ki2 

0 


^12 

-(K12+K23) 


K 


23 


0 

^23 

-(Kgg+Kg^) 


K 


34 


K 


34 


-K 


34 


T = 


and 


r 





F = 


where, 


P 

V 

H 

h 

A 



r H T 


AW 





= mass density 

= volume associated with node 
= specific heat- 

= convective coupling, hA 
= convective heat transfer coefficient 
= heat transfer area 

= conductive coupling between nodes i and 3 , kA/AX.j 


k = thermal conductivity 

AX. . = conduction distance between nodes i and j 

13 

T adiabatic wall temperature of fluid 

AW 
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To obtain the solution to theEq. (F-1) , we premultiply both sides of (F-1) by 
the inverse of M and obtain 

T = m“^T + M'^F Eq. (F-2) 

If we define a symmetric matrix A by 

A = Eq. (F-3) 

— 1/2 1/2 

we obtain after premultiplication of Eq. (F-3) by M*" ' and postmultiplication by M ' 

Eq. (F-4) 

After substituting Eq. (F-4) into Eq. (F-2) and defining a matrix N for convenience in 
writing by 


N = Eq. (F-5) 

we get 

T = NT + M~^F Eq. (F-6) 

We now assume that the solution toEq. (F-1) can be expressed as 

T= 0 (t) + Too Eq. (F-7) 

where 0 is a vector having time-dependent components andT^o is a vector that is 
independent of time. (Physically, T represents the steady-state temperatures that 

— CO — 

the system will achieve. ) 

Substitution ofEq. (F-7) into Eq. (F-6) yields; 

0 = N0 +NT + M“^F Eq. (F-8) 

— ^ ^00 
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Eq. (F-9) 


and we shall let T be such that It satisfies 


NT +M“^F=0 
—06 — 


Therefore, we have 


•^CO — 


Then with satisfying Eq. (F-9), Eq. (F-8) becomes 


e =Ne 


which has the solution*’'' 


e = e*N 

— -o 


Returning to Eq. (F-7) , we have after substituting Eq. (F-11) 


T = e^^ 9 + T 

— - o — < 


To evaluate 9^ we use the initial condition 


T(0) = Tj 


thus 


e =: T, - T 
— o — i — oo 


which gives upon substitution into eq, (F-12) 

T = Tt. - T 1 + T 


Eq. (F-10) 


Eq. (F-11) 


Eq, (F-12) 


Eq. (F-13) 


*See for example Hochstadt, "Differential Equations," Holt, Rinehart and Winston, 
pp 55-58. 
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Eq. (F-13) is the solution to Eq, (F-1) along with the Initial condition; 
however, from a computational standpoint Eq. (F-13) is not convenient because 
evaluation of the e^qponential term requires summing an infinite power series in 
the matrix N, 

We now develop a more convenient form for Eq. (F-13) . Recalling the definition 
of N from Eq. (F-5) 


N = M 


- 1/2 


AM 


1/2 


we see that 


and, in general, that 


Therefore, the Taylor series 
00 

tN_ 

" "2 

J=o 

can be written 

etN = M-l/2 I ^ j^jl/2 
j=0 

which leads to the result 

_fN _ _ „-l/2 tA 

e — e = M e M 



Eq. (F-14) 


Because the matrix A is a symmetric matrix we can find a matrix V such that 


V“^ AV =A 


Eq. (F-15) 
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where A Is a diagonal matrix formed with the eigenvalues of A? The matrix V can be 
selected to be orthogonal, that Is, It Is a matrix whose columns comprise the elements 
of n linearly independent eigenvectors of matrix A that are mutually orthogonal and of 
unit length. 

From Eq. (F-15) we have 
A = V A 


which when substituted into Eq. (F-14) yields 


tN _ ,,-1/2 tA 
e = M V e 


V-l 


Eq. (F-16) 


where we have made use of the identity 

tVAV^ _ tA 
e = Ve V 

tN 

The representation for e given in Eq. (F-16) is convenient since A is a 
diagonal matrix and therefore the exponential term is easily and explicitly evaluated 
as follows: 


tA _ 


e = 




tA2 



0 


e^^n 


where the X's are the eigenvalues of matrix A. 


*cf. Hildebrand, "Methods of Applied Mathematics Prentice-hall , pp. 37-39. 
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Substituting Eq. (F-16) into Eq. (F-13) we obtain the final form for the 
solution to equation (F-1) 


T = r^i “ T 1 + T Eq. (F-17) 

where we have made use of the property that for normalized modal matrices , the trans- 
pose of the matrix equals the inverse of the matrix, i.e. , 



Equation (F-17) is Eq. (9) of NASA CR-2435 by Maise and Rossi; it forms the 
basis of the HAN method. 

It should be noted that for a system with n nodes or elements the matrix V will 
have n columns of eigenvectors and the matrix A will have n eigenvalues alor^ the 
diagonal. As noted in NASA CR-2435, good approximate solutions for T are obtained 
by finding and using only the "dominant" eigenvectors and eigenvalues of A. Large 
savings in computer time (factor of ten or more) are achieved by finding and using 
only say 5 dominant eigenvectors and eigenvalues for a system of 100 nodes or more. 
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